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CHAPTER  1 


INTRODUCTION 


Predicting  the  chemical  composition  of  a multi-phase,  multi- 
component  system  at  equilibrium  is  of  much  practical,  as  well  as 
theoretical,  interest.  This  can  be  posed  as  a nonlinear  programming 
problem  in  which  a convex  function,  called  the  free  energy,  is  minimized, 
subject  to  linear  mass  balance  equations  in  nonnegative  variables. 

There  is  an  associated  dual  problem,  equivalent  to  a geometric  program. 
The  main  purpose  of  this  work  is  to  study  this  "chemical  duality"  and 
apply  the  existing  theory  of  geometric  programming  to  analyze  and  solve 
chemical  problems . 

Willard  Gibbs  £30],  who  developed  in  1876  the  fundamental  concept 
of  free  energy,  paved  the  way  for  a systematic  study  of  chemical 
equilibrium.  For  many  years  the  inability  to  handle  large  complex 
systems  manually  and  lack  of  theory  dealing  with  minimization  problems, 
prevented  the  use  of  his  powerful  elegant  ideas. 

In  the  last  thirty  years,  rapid  developments  of  both  optimization 
theory  and  computers  proved  the  usefulness  of  Gibbs'  ideas  and  their 


superiority  over  earlier  methods  involving  solution  of  nonlinear 


L 


equations.  Following  the  work  of  Brinkley  [10],  a serious  effort 
began  at  the  RAND  Corporation  in  the  fifties,  contributing  applica- 
tions, theory,  and  computational  methods  Notable  among  these  are 
the  works  of  Bigelow  [8,9],  Clasen  [12-14],  Dantzig  [l6-l8],  and 
Shapiro  [44-47],  whose  theoretical  developments  form  a framework 
to  our  work. 

Geometric  programming  theory  was  developed  in  the  sixties  by 
Duffin,  Peterson  and  Zener  [25] . The  relation  of  Geometric  Programming 
to  chemical  equilibrium  was  shown  by  Avriel  [4],  and  Passy  and  Wilde 
[40-41] . 

In  addition  to  well  known  applications  of  chemical  equilibrium 
to  the  study  of  chemical  reactions  in  chemical  process  design,  rocket 
propellants,  and  the  oil  industry,  there  are  potentially  many  applica- 
tions in  biochemistry  and  medicine.  Biological  systems  are  in  many 
cases  closer  to  chemical  "ideality"  than  are  commercial  inorganic 
processes,  but  their  size  and  complexity  were  until  recently  beyond 
solution  by  available  computational  and  analytical  methods.  Very 
little  is  known  today  about  the  intricate  chemistry  of  many  physi- 
ological functions,  for  example,  the  way  medications  react  in  the  human 
body.  This  work  is  motivated  in  part  by  the  need  to  analyze  these 
and  related  systems. 

Chapter  2 reviews  the  formulation  of  the  problem.  Chapter  5 
studies  thermodynamic  concepts,  deriving  general  properties  of  free 
energy  functions  from  basic  principles.  Geometric  programming  and 


/ 
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the  dual  chemical  problem  aie  introduced  in  Chapter  4,  with  a new 
chemical  interpretation  of  the  dual  variables.  In  Chapter  5 we  apply 
the  chemical  duality  and  the  results  of  Chapter  4 to  study  the  nature 
of  equilibrium  solutions.  It  is  shown  that  uniqueness  of  solutions 
is  insured  under  regularity  conditions  relating  to  the  dual  problem. 
Geometric  programming  methods  can  also  handle  degenerate  systems  and 
degenerate  solutions . 

The  theoretical  discussion  is  followed,  in  Chapter  6,  by 
several  applications  of  duality  in  characterization  and  formulation 
of  chemical  problems.  Among  the  results  is  a generalization  of  a 
"Goaling"  method,  in  which  side  conditions  are  added  to  the  problem 
and  some  of  its  mass  balance  equations  are  relaxed. 

Chapter  7 presents  a dual  algorithm  for  solving  the  chemical 
equilibrium  problem.  The  method,  adapted  from  an  algorithm  proposed  by 
Dembo  [20],  is  shown  to  be  a convex  cutting  plane  algorithm.  It  solves 
a transformed  geometric  program  through  a sequence  of  "condensed" 
linearized  programs.  Comparative  results  of  tests  with  this  algorithm 
are  presented,  together  with  a small  collection  of  test  problems, 
which  appear  in  an  appendix. 

The  last  chapter  presents  a new  extension  to  geometric  pro- 
gramming, which  includes  variables  as  exponents.  This  extension, 
called  "transcendental  geometric  programming  was  primarily  developed 
because  of  its  potential  applications  to  engineering  problems.  How- 
ever, a dual  problem  to  the  transcendental  program  is  shown  to  be  an 
interesting  generalization  of  the  chemical  equilibrium  problem,  one 
particularly  suitable  for  treating  nonideal  systems. 
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CHAPTER  2 


FORMULATION  AND  BASIC  CONCEPTS 

This  chapter  introduces  notation  and  terminology  used  through- 
out this  work.  The  chemical  system  is  described,  together  with  the 
concepts  required  for  mathematical  formulation  of  the  chemical 
equilibrium  problem. 

2.1  General  Notation  Conventions 

We  denote  by  1 the  real  line;  En  an  n-dimensional  Euclidean 
space,  and  E^  its  positive  orthant.  There  will  be  no  distinction 
in  notation  between  scalars  and  vectors  since  the  dimensions  will  be 
clear  from  the  context.  Vectors  are  assumed  to  be  column  vectors. 

The  transpose  of  a vector  x is  denoted  x'.  However,  the  prime  will 
be  omitted  at  times,  when  no  ambiguity  arises.  Inner  products  will  be 
written  xy  and  sometimes  x' -y,  where  x and  y are  vectors.  Matrices 
will  always  be  denoted  by  capitals.  If  A is  a matrix,  A,  will 
denote  its  j-th  column;  A.  its  i-th  row.  The  symbol  (k)  denotes 
a set  of  consecutive  integers  belonging  to  the  k-th  partition  of  some 
integer  set  N. 


k 


References  are  enclosed  in  square  brackets  [ ] . The  material 
is  divided  into  chapters;  the  chapters  into  sections.  Theorems,  lemmas, 
propositions  and  other  major  items  are  numbered  sequentially  within  ea_h 
chapter  in  the  form  k.n  where  k is  the  chapter  number.  Other 
items  and  equations  are  numbered  sequentially  within  each  section  by 
a single  number.  Internal  reference  to  equation  i is  by  (i)  if 
the  equation  is  in  the  same  section,  (j.i)  if  it  is  in  section  j of 
the  same  chapter  and  (k.j.i)  if  it  is  in  chapter  k section  j. 


2.2  The  Chemical  System 

The  description  below  is  based  mostly  on  the  notation  of 
Shapiro  and  Shapley  [47] . Some  changes  were  introduced  to  adapt  the 
formulation  to  Geometric  Programming  conventions. 

The  system  under  consideration  is  composed  of  a finite  number 
(K)  of  homogeneous  phases  denoted  4^,  4^,  ...  , 4^.  We  shall  some- 
times refer  to  phase  4^  simply  as  phase  k.  A phase  is  homogeneous 
when  every  part  of  it  has  the  same  pressure,  temperature  and  chemical 
composition.  By  this  definition,  a homogeneous  phase  need  not 
occupy  a contiguous  space  in  the  system.  For  example,  oil  drops  in 
water  can  be  considered  a homogeneous  phase. 

Each  phase  is  a subset  of  a set  of  chemical  species . The 


species  are  those  chemical  compounds  which  one  can  expect  to  find  in 
in  the  system.  The  question  of  which  compounds  should  be  considered 


when  constructing  a model  is  best  answered  on  the  basis  of  experience  and 

judgment.  In  general  the  choice  of  species  is  up  to  the  designer 

of  the  model.  The  mathematical  analysis  of  the  system  assumes  a 

given  set  of  species  S^,  Sg,  ...  , Sn,  and  the  results  may  be 

mathematically  valid  even  if  the  system  is  chemically  unsound. 

For  our  purposes  the  species  will  be  partitioned  in  the  phases, 

that  is,  each  species  belongs  to  one,  and  only  one  phase.  Thus  a 

chemical  compound  having  the  same  chemical  formula  but  appearing  in 

two  different  phases  is  considered  as  two  different  species.  We 

shall  denote  the  fact  that  species  S.  appears  in  phase  ® by 

J & 

S 6 t,  or  equivalently  j £ (k),  where 
J k 


(k)  = tj|Sj  e *k) 


denotes  the  set  of  indices  of  species  in  . The  phase  to  which 

Sj  belongs"  is  denoted  by  ^us  we  can  write  k(j)  = {k|  j £ (k)}. 

If  i £ (k( j))  then  S.  and  S belong  to  the  same  phase. 

1 J 

The  quantity  of  species  S.  in  the  system  is  denoted  by  x , 

J J 

measured  in  moles . (x,  is  thus  proportional  to  the  number  of  molecules 

d 

of  S,  in  the  system.)  The  values  of  x,  are  usually  unknown  and  part 
J J 

of  the  problem  to  be  solved  is  to  find  the  vector  x = '{x^,  Xg,  ...  , xr) 
of  chemical  composition  when  the  system  is  at  equilibrium.  The  com- 
position-vector x Is  partitioned  into  the  phases  in  the  same  way  as 
described  for  the  corresponding  species,  so  that  we  can  describe  the 


6 


For  simplicity  we  will 


composition  of  phase  by  the  vector  x^. 
omit  the  brackets  and  write  x^.  It  will  always  be  clear  from  the  con- 
text whether  the  reference  is  to  a vector  or  to  a single  component  of 
a vector. 

We  define  the  integer  sets 

(1,  2,  ...  , K}  (the  phase  indices) 

N = {1,  2,  ...  , n)  (the  species  indices) 

The  total  moles  of  all  species  in  a phase  is  denoted  by 


x,  = £ x. 

k j€<k>  J 


k 6 K 


(1) 


and  the  mole  fraction  of  species  (in  is  denoted  by 


WVj)  "hen  “d  *k(j) 


> 0 


(2) 


The  mole  fraction  serves  as  a dimensionless  concentration.  Note  that 
the  concentration  x^  is  defined  only  when  x^^  > 0,  namely,  when 
the  phase  actually  exists. 

Each  of  the  species  in  the  system  is  composed  of  a set  of 


basic  units  called  subspecies  which  are  usually  smaller  chemical  units 
like  atoms,  ions  or  radicals.  In  general,  the  subspecies  are  units 


which  do  not  decompose  into  any  smaller  units  in  the  given  system. 


The  choice  of  subspecies,  much  like  the  selection  of  species  to 

be  present  in  the  model,  is  a matter  of  experience  and  judgment. 

However,  once  the  species  are  specified,  the  set  of  subspecies  must 

provide  for  every  species  so  that  each  species  is  composed  of  one 

or  more  of  the  given  subspecies,  denoted  by  B^,  Eg,  ...  , Bm> 

The  composition  of  one  mole  of  species  S.  can  now  be  described 

J 

by  the  "equation" 


m 

£ Vi  - sj  J £ n (5) 


here  a.  is  the  number  of  moles  of  B.  in  one  mole  of  S . The 

J-  J J 

constants  a.  form  a matrix  A 6 ^nXn.  The  elements  of  A need 
J 

not  all  be  nonnegative,  although  in  simple  cases  they  usually  are. 


With  each  species  S . we  can  associate  a column  vector  of  A as 

J 

indicated  by  (j).  Hus  vector  A is  called  the  formula  vector 

J 

of  S . 

J 


With  the  notation  developed  so  far  we  can  now  describe  the 
chemical  structure  of  the  system  by  the  triplet  (B,S,4>)  where  B 
is  an  m-vector  of  subspecies,  S is  an  n -vector  of  species  and  $ 
is  a K- vector  of  phases  which  induces  a partition  of  the  set  N. 
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2.3  The  Mass  Balance  Constraints 


Except  where  noted  otherwise,  we  shall  assume  that  the  chemical 
system  is  closed,  that  is,  there  is  no  flow  of  matter  in  or  out  of  the 
system.  The  contents  of  the  system  are  thus  determined  by  the  initial 
input  to  it,  usually  measured  in  terms  of  moles  b^  of  the  subspecies 

V 

Although  the  actual  input  is  usually  implemented  by  introducing 
certain  species,  nothing  is  lost  by  assuming  that  these  were  decomposed 
to  their  subspecies  which  were  then  separately  introduced.  Species 
can  decompose  or  form  via  chemical  reactions  and  there  is  no  law  of 
conservation  for  moles  of  species.  Subspecies  on  the  other  hand,  by 
their  definition  cannot  form  or  decompose  so  their  quantities  must 
be  conserved,  regardless  of  the  species  in  which  they  appear. 

Given  the  initial  amount  ^ of  B^,  i € M = (1,  2,  . . . , m) 
we  can  write  for  any  composition  vector  x 

n 

Z a.  x.  = b.  i £ M (l) 

j=l 

These  relations,  written  in  matrix  form 

Ax  = b 

are  called  the  mass  balance  constraints . The  composition  vector  x 
cannot  have  any  negative  components  so  we  require  naturally  that 
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x > 0 


Some  remarks  concerning  the  matrix  A are  in  order  here. 
Although  a^j  is  the  "amount  of  in  " it  need  not  represent 

an  actual  physical  quantity.  It  reflects  only  the  structure  of  S 

J 

with  the  given  set  of  B^’s.  Sometimes  it  is  convenient  to  define 

the  B 's  in  a way  which  assigns  a negative  value  to  a for  some 
I JL  J 

l and  j.  In  these  situations,  the  rest  of  the  B. 's  are  defined 
in  such  a way  that  B^  is  implicitly  contained  in  some  of  them,  so 
that  equation  (3)  is  more  of  an  accounting  device  than  a meaningful 
physical  relation.  In  some  models  additional  constraints  reflecting 


electroneutrality  are  included  in  the  matrix  A.  This  has  the  effect 
of  introducing  a new  subspecies  (electrical  charge)  which  normally 
has  a quantity  of  0 on  input  (i.e.,  b^  = 0 for  the  charge). 

The  previous  remarks  notwithstanding,  one  can  always  generate 
completely  arbitrary  matrices  A,  and  associate  with  them  abstract 
chemical  systems , as  opposed  to  real  chemical  systems ■ We  shall  be 
interested  mostly  in  the  latter,  but  the  former  have  useful  appli- 


cations since  they  are  dual  to  geometric  programs,  as  will  be  shown 
in  Chapter  h. 
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2.k  Free  Energy  Functions  and  Ideality 


Although  most  of  our  chemical  and  mathematical  ingredients  are 
at  hand,  nothing  has  been  said  so  far  about  the  mathematical  formula- 
tion and  characterization  of  equilibrium.  Traditionally,  equilibrium 
conditions  were  defined  in  terms  of  equilibrium  constants  and  the 
"Mass  Action  Laws"  to  be  described  in  the  following  section.  Gibbs 
[30]  defined  a potential  function  F related  to  the  chemical 
potential  of  a system  and  showed  that  equilibrium  is  attained  when 
F reaches  its  minimum.  The  function  F,  called  the  free  energy 
function, will  be  described  in  more  detail  in  the  next  chapter.  For 
the  time  being  we  note  that  in  general,  F has  the  form 


Here 


n 

F(T,P,x)  =RT  Z x n (T,P,x) 
j=l  J J 


(1) 


P 

R 

T 

x 

Hj(T,P,x) 


is  the  total  pressure 
is  the  gas  constant 
is  the  absolute  temperature 
is  the  composition  vector 

is  a function  called  the  chemical  potential  of 


For  most  of  our  analysis  we  assume  that  the  system  is  maintained 
under  fixed  pressure  and  temperature,  so  that  T and  P are  constants. 
For  convenience  we  then  eliminate  the  constant  RT  and  redefine  the 
dimensionless  free  energy  at  fixed  T and  P 
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(2) 


n 

F(x)  = £ x Li  (x) 
j=l  J J 

The  assumption  of  fixed  temperature  and  pressure  is  usually  made 
when  the  system  is  free  to  expand  into  the  atmosphere 
surrounding  it  and  is  in  contact  with  a large  heat  reservoir 
(e.g.  the  surrounding  air).  This  is  especially  true  of  biological 
systems  which  maintain  constant  temperatures  through  internal 
mechanisms,  and  cannot  tolerate  large  variations  in  either  temper- 
ature or  pressure. 

From  the  form  (2)  of  the  free  energy  function  it  is  clear 
that  it  can  be  written  as  the  sum  of  free  energies  of  the  species 
in  the  system.  Denoting  the  free  energy  of  Sj  by  f^(x)  we  write 


F(x)  = 


E f.(x) 

j=i  J 


and 


fjW 


Vj(x) 


The  function  n (x)  can  now  be  interpreted  as  the  partial  molar 

free  energy  of  S • For  x > 0 we  can  assume  with  no  loss  of  generality 
J J 

that  for  any  P,  T and  x 


Uj(P,T,x)  = c j(P,T,x)  + log  Xj 


(3) 


I 

\ 
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The  functions  Cj(P,T,x)  can  (and  often  do)  have  a rather  complex 

form.  However,  much  of  our  attention  will  be  focused  on  ideal 

systems  in  which  Cj(P,T,x)  is  independent  of  x.  Under  the  hypothesis 

of  fixed  P and  T,  the  functions  c.(*)  become  constants  c , and 

J J 

we  have  for  x > 0 


F(x)  = E xj(cj  + lo6  xj) 

J 1 


(V 


The  assumption  of  ideality  is  by  no  means  a mere  theoretical 
construction.  Many  systems,  especially  those  under  fairly  low 
pressures  and  moderate  temperatures,  exhibit  ideal  behavior  over  some 
”ange  of  P,  T and  x.  Most  biological  systems  satisfy  these  conditions 
and  can  therefore  be  treated  as  ideal  systems.  Dantzig  et  al.  [17], 
who  modelled  the  human  respiratory  system,  obtained  results  which 
are  in  excellent  agreement  with  experimental  data  by  using  an  ideal 
free  energy  of  the  form  (4). 

Equation  (3)  is  not  the  only  form  expressing  the  partial  molar 
free  energy  [44] . In  chemical  thermodynamics  the  most  widely  used 
(dimensional)  equivalent  of  (j)  is 


hj(P,T,x)  = dj(P,T)  + RT  log  aj 


(5) 


where  a,  is  the  activity  of  S , which  usually  depends  on  P,  T 
J J 

and  x.  The  bar  indicates  that  p is  a dimensional  function. 
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Dividing  (5)  by  RT,  it  is  easy  to  verify  that  our  previous  definition 
of  ideality  implies 


*j  = r (p,*)-«j 


(6) 


Here  y,  is  called  the  activity  coefficient.  In  the  general 
J 

(nonideal)  case  y,  is  a function  of  P,  T and  x.  Some  of  the 
J 

more  fundamental  properties  of  free  energy  functions  will  be  developed 
and  analyzed  in  the  next  chapter. 


2.5  The  Chemical  Equilibrium  Problem 

We  can  now  state  the  mathematical  form  of  the  chemical 
equilibrium  problem,  assuming  that  the  system  is  ideal  and  under 
fixed  pressure  and  temperature. 

Problem  CPI 

Minimize  F(x)  = x‘*(  + log  x)  (l) 

subject  to  Ax  = b 
x > 0 

Here  x £ in,  c £ £n,  A £ t?1*11 , b £ M*'.  The  notation  log  x represents 
a vector  with  components  log  Xj • When  x^  = 0,  log  x^  is  undefined. 

To  maintain  continuity  of  F over  the  nonnegative  orthant  of  B11, 
we  define 
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x log  x = 0 when  x = 0 . 
J J J 


With  this  definition  F is  well  defined  and ''right  continuous"  over 
its  domain.  Note  that  £ implicitly  partitions  the  species  into 
the  phases.  The  matrix  A is  partitioned  into  submatrices 


A(k)  = (Aj  I J £ 


k € K 


A^)  contains  the  formula  vectors  of  species  in  phase  4>^. 

In  the  general  case,  under  fixed  temperature  and  pressure  the 
problem  is 


Problem  CPN 


minimize  F(x)  = x'  n(x) 

subject  to  Ax  = b 


x > 0 


where  we  define  = ® when  x^  = 0 to  maintain  continuity. 

It  is  assumed  that  p (x)  is  welJ  defined  for  x > 0. 

J J 

The  problem  of  finding  the  equilibrium  composition  for  a system 
{B,S,$>)  can  now  be  represented  by  the  triplet  {F,A,b)  where  F is 
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the  free  energy  function,  A is  a (partitioned)  m X n matrix  and  b 


M/S 


is  an  m- vector. 


We  define  the  set  of  possible  solutions 


X * X(A,b)  2 [x  € ^(Ax  = b,  x > 0) 


A composition -vector  x is  said  to  be  feasible  if  x € X. 


■uwyiww,1"?'"-  J"  ■ 


CHAPTER  J 


THERMODYNAMICS  OF  CHEMICAL  EQUILIBRIA 


There  are  two  basic  approaches  to  the  study  of  chemical 
equilibrium;  one  dealing  with  processes,  the  other  dealing  with 
states . In  this  chapter  we  shall  first  review  some  notions  of 
the  classical  "process"  approach  and  then  develop  general  properties 
of  free  energy  functions  based  on  the  "state"  approach. 


3.1  Reactions  and  Reaction  Vectors 

Chemical  species  usually  react  to  produce  other  species. 

In  principle,  all  chemical  reactions  taking  place  in  a system  proceed 
simultaneously,  each  at  some  characteristic  reaction  rate  which 
depends  on  temperature,  pressure  and  concentrations  of  the  species 
in  the  system.  Equilibrium  is  reached  when  no  observable  change 
with  time  occurs  in  the  system,  in  other  words,  when  reaction  rates 
are  all  zero. 
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A typical  chemical  reaction  can  be  expressed  schematically  by 


E r?j  <_ 

j £ R J J 


l PlS 


j£P 


J J 


(i) 


Here  R and  P are  the  sets  of  reactants  and  products  respectively; 


r'  and  p'  are  the  stoichiometric  coefficients  of  the  reactants 
J J 


and  products  respectively  The  double  arrow  indicates  a possible 
reaction  in  the  reverse  direction. 

Letting 


'J 


rJ  = 


0 


j £ r 

j 4 R 


(2) 


.h 


’ 0 


J 6 P 

J 4 p 


(5) 


we  can  write 


n * n 

t rS  < E p S 

j=l  J 0 J=1  J 


j 


00 


The  relation  above  can  be  expressed  as  a formal  equation  by  noting 

that  Sj  is  represented  by  a column  vector  A.  via  (2.2.3).  Let 

J 


6 = p - r,  for  all  j £ N.  From  (4)  we  get 

J J J 


1 
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E 0.S.  < — o 

j=i  J J 


which  can  be  converted  to  equations  of  mass  conservation 


E 

J=i 


= o 


or  for  short 


A0  = 0 


(5) 


The  vectors  0 £ ffin  satisfying  (5)  are  called  reaction  vectors . 

By  this  definition,  the  set  of  all  possible  reaction  vectors  is  the 
null  space  of  A. 

Not  all  possible  reactions  implied  by  (5)  can  occur  in  a 


system  with  a given  composition  x.  If  x = 0,  species  S cannot 

J J 


be  a reactant,  so  that  only  reaction  vectors  having  0 > 0 are 

J 


possible.  We  are  thus  led  to  the  definition  of  feasible  reaction 
vec  ors  at  a given  composition  x,  as  those  reaction  vectors  0, 


for  which  0,  > 0 whenever  x = 0. 

J J 


For  x 6 X and  0 feasible  at  x we  can  find  e > 0 
sufficiently  small  so  that  for  | 0 | < e 

A(x  +0)  = b,  x + 0 > 0 

Hence 


x + 0 C X 


This  observation  identifies  feasible  reaction  vectors  as  feasible 


directions  for  changes  in  composition.  When  a reaction  characterized 
by  0 takes  place,  x changes  in  the  direction  9.  A system  for 
which  there  exists  at  least  one  nonzero  vector  9 satisfying  (5) 
is  said  to  be  reactive . This  work  deals  only  with  reactive 
systems . 

As  mentioned  earlier,  there  is  no  loss  of  generality  in 
assuming  that  the  subspecies  II  are  elements,  in  which  case  A > 0. 
(For  the  moment  we  ignore  electroneutrality  equations.)  With  this 
assumption  we  can  state  the  following  intuitive  result. 

Proposition  5.1:  For  a real  chemical  system,  the  set 

X s (x  € Kn|Ax  = b,  x > 0}  is  bounded. 

Proof:  We  may  assume  that  X is  nonempty,  and  that  A > 0. 

If  X were  unbounded,  there  would  exist  an  x^  £ X and  a 
nonzero  reaction  vector  9 > 0 such  that 

x^  + 7\9  <Z  X for  all  X > 0 

Since  the  matrix  A would  have  at  least  one  positive  element  in  each 
column,  the  equation  A0  - 0 would  imply  0 = 0,  a contradiction.  □ 

This  result  will  be  useful  in  tne  study  of  properties  of 
equilibrium  sets. 
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Corollary  5.2:  Any  nonzero  reaction  vector  must  have  at  least  one 

positive  and  one  negative  component. 


Intuitively,  the  corollay  states  that  a species  can  form  by 
reaction  only  if  at  least  one  other  species  decomposes. 


3-2  The  Mass  Action  Lavs 

The  classical  approach  to  chemical  equilibrium  considers 
reactions  as  "reversible"  processes  proceeding  both  forward  and  back- 
ward simultaneously  [31] . For  the  reaction 


z rj3j  ■ z pjsj 

j=i  J J j=i  J J 


the  forward  rate  pf  is  given  by 


n r , 


Pf  = k (r)  n x * 
j=l  J 


Similarly,  the  backward  rate  is 


n P. 

pb  = Mp)  n 

0=1  J 


Here  k^(r)  and  k^(p)  are  the  rate  constants  for  the  forward  and 
backward  reactions  respectively,  the  system  being  assumed  ideal. 
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Equilibrium  is  attained  when  the  forward  and  backward  rates 
are  equal,  so  that  no  net  change  is  observed  with  time.  Therefore, 
at  equilibrium  we  must  have 


k1(r)  n x.J  = k2(p)  II 


Using  the  definition  of  8 from  the  previous  section  we  have  at 
equilibrium 


n 

n 


J=1 


k2(p) 

= k^rj 


H K(0 ) 


(4) 


Equation  (4)  is  called  the  mass  action  law.  It  has  to  be 
satisfied  for  every  feasible  reaction  vector  6 at  equilibrium.  It 
can  be  shown  [7,47]  that  (4)  is  a necessary,  as  well  as  sufficient 
condition  for  x to  be  an  equilibrium  composition  of  an  ideal 
system  when  x £ X. 

The  mass  action  laws  are  closely  related  to  the  Kuhn-Tucker 
Conditions  [36]  for  the  minimum  of  F(x)  in  problem  CPI.  We  shall 
show  in  the  next  chapter  that  (4)  appears  in  a geometric  program 
associated  with  CPI. 

When  the  matrix  A has  full  rank — (m),  the  dimension  of  its 
null  space  is  n-m.  Thus  there  are  no  more  than  n-m  linearly  inde- 
pendent reaction  vectors . 

Taking  the  logarithm  of  (4)  we  have 
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(5) 


9 log  x = log  K(0) 

\ 

where  0j  log  x^  = 0 when  x^  = 0. 

Equation  (5)  indicates  that  log  K (0)  is  a linear  function 
of  9,  and  therefore  (5)  represents  a linear  system  of  equations 
in  x with  coefficients  9 from  the  null  space  of  A.  By  the  pre- 
ceding remark,  this  system  has  at  most  n-m  independent  equations, 
which  together  with  the  mass  balance  equations  (2.3.I)  determine  x. 

3.3  Systems,  Properties,  and  States 

Much  of  the  terminology  used  in  thermodynamics  refers  to 
observable  physical  entities  and  to  measurable  quantities.  Thus  many 
of  the  definitions  are  operational  rather  than  conceptual. 

Hie  following  definitions,  though  not  completely  precise 
mathematically,  will  serve  our  purposes. 

(i)  A system  is  any  collection  of  well  defined  physical  entities, 

(ii)  A property  is  any  measure  of  the  system. 

(iii)  A state  of  a system  is  the  collection  of  values  of  all  its 
properties. 

For  a given  chemical  system  the  temperature  and  pressure  are 
properties,  as  are  the  location,  direction,  and  velocity  of  every 
molecule  in  the  system.  It  would  seem  that  a specification  of  the 
state  would  require  an  almost  infinite  amount  of  information. 


23 


Fortunately,  only  a limited  number  of  properties  are  of  interest  in 
thermodynamic  studies.  Moreover,  experience  has  shown  that  these 
properties  are  not  independent  of  each  other.  For  most  purposes 
then,  only  a small  set  of  independent  properties  need  be  specified 
in  order  to  determine  all  other  properties  of  interest  and  the  state. 
These  observations  lead  to  formulating  equations  of  state  which 
relate  independent  and  dependent  properties.  In  what  follows  we 
occasionally  use  the  term  variables  for  properties.  There  are 
two  types  of  variables: 

(i)  Intensive  variables  (e.g.  pressure  P and  temperature  T) 
which  are  independent  of  the  mass  of  the  system. 

(ii)  Extensive  variables  (e.g.  volume  V,  energy  U,  composition  x) 
which  are  mass-dependent.  It  will  be  shown  later  that  this 
classification  is  significant  in  duality  relations. 

One  point  to  be  emphasized  about  properties  and  states  is 
that  they  refer  to  points  stationary  in  time.  It  is  therefore 
appropriate  to  speak  of  properties  as  "state  variables, " in  contrast 
to  quantities  measured  during  transition  like  heat  flow  into  the 
system  or  work  done  by  the  system. 


3.4  Spontaneous  Processes  and  the  Gibbs  Free  Energy 

In  the  study  of  chemical  processes,  one  often  faces  the  question: 
"Will  this  process  (or  reaction)  occur  spontaneously  under  given 
conditions?".  The  second  law  of  thermodynamics  provides  a quantitative 
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measure,  the  entropy  function  S,  to  answer  this  question.  One  con- 


sequence of  this  law  is  that  every  spontaneous  process  is  accompanied 
by  a net  increase  in  entropy  for  the  system  and  its  surroundings. 

This  criterion,  though  theoretically  valid,  is  nevertheless  of  little 
practical  use  since  it  requires  some  knowledge  of  the  surroundings. 

The  free  energy  functions  F and  A introduced  by  Gibbs  and  Helmholz 
resolve  this  difficulty.  In  terms  of  other  thermodynamic  properties, 

F is  defined  by 


F - U + PV  - TS 


(1) 


Here  U is  the  internal  energy  of  the  system,  P,  T,  and  V are 
the  pressure,  temperature,  and  volume  respectively.  Similarly, 


A = U - TS 


<2) 


Note  that  both  F and  A are  extensive  variables  since  U and  S are. 

We  shall  assume  that  the  system  is  closed,  under  fixed  temper- 
ature T and  pressure  P,  and  that  it  undergoes  a process  from  state  1 
(Ul,  Vi,  Si)  to  state  2 (U2,  V2>  S2).  Following  Denbigh  [21],  we 

let  q be  the  heat  absorbed  by  the  system;  w,  the  work  done  by  it  in 
the  process.  The  first  law  of  thermodynamics  states  that 


U2  " U1  = q " W 


(3) 


The  second  law  requires  that  for  every  spontaneous  process 


q < T(S2  - Sx) 
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From  ( 1 ) , ( 3 ) we  have 


Fg  - Fx  - q - w + P(V2  - V]_)  - T(S2  - S1) 


And  by  (4) 


F2  - Fx  < - w + P(V9  - Vn) 


2 


The  work  w can  be  partitioned  into  pressure -volume  work  P(v^  - V ) 
and  other  (chemical)  work  w1  so  that 


w = P(V2  - Vx)  + w' 


p _ p < - y 1 
r 2 1 - 


Since  no  external  work  is  done  on  the  system,  w 1 > 0,  whence 


F,  - Fx  < 0 


The  inequality  in  (4)  and  (6)  is  strict  when  the  process  is 
irreversible  as  are  all  real  processes.  We  conclude,  therefore,  that 
every  spontaneous  process  in  a closed  system  at  fixed  T and  P is 
accompanied  by  a decrease  in  its  free  energy.  (A  similar  analysis 
holds  for  A under  fixed  T and  V.) 

The  free  energy  during  the  process  can  be  described  as  a 
function  of  time  t from  the  start  of  the  process.  F(t)  is  thus 
a decreasing  function.  Taking  the  limit 


lim 

At  -»  0 


F(t  + At)  - F(t)  _ dF(t) 
’ At  " dt 
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we  conclude  from  (6)  that  for  a real,  spontaneous  process 


0 


(7) 


This  conclusion  is  extensively  treated  by  Bigelow  [7],  who  calls  it 
the  "global  least  action  principle. " 


3.5  Equilibrium  States 

From  a theoretical,  as  well  as  computational  point  of  view, 
there  are  several  advantages  to  the  "state"  approach  over  the  "process" 

• '! 

{ approach.  Following  Callen  [11],  we  present  here  such  a state  approach 

and  derive  some  new  mathematical  consequences  of  its  postulates. 

Hie  first  and  second  laws  of  thermodynamics  will  be  replaced 
here  by  two  postulates  on  simple  systems . A simple  system  is  a homo- 
geneous, uncharged,  and  chemically- inert  system  which  is  not  acted 
upon  by  any  external  fields  like  gravity  or  magnetic  fields . 

Postulate  1:  There  exist  states,  called  equilibrium  states,  of 

simple  systems  in  which  the  system  is  completely  characterized 
macroscopically  by  its  internal  energy  U,  its  volume  V and  its 
composition  vector  x. 

For  the  second  postulate  we  define  a composite  system  to  be 
a collection  of  two  or  more  simple  systems . A composite  system  is 
said  to  be  closed  if  no  flow  of  energy  or  matter  can  occur  through 
its  boundaries  and  if  its  volume  is  fixed.  All  restrictions  on 

I 
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flow  of  energy  or  matter  between  the  simple  systems  in  a closed 
composite  system  are  called  internal  constraints. 

Postulate  2:  There  exists  a function  S = S(U,V,x),  called  the  entropy 

of  a composite  system,  defined  for  all  equilibrium  states  and  having 
the  following  properties: 

(i)  S is  differentiable 

(ii)  S is  monotone  increasing  in  U 

(iii)  S is  additive  over  the  constituent  simple  systems 

(iv)  The  values  assumed  by  the  extensive  variables  (U,V,x) 

in  the  absence  of  internal  constraints,  maximize  S over  the 
manifold  of  internally  constrained  equilibrium  states. 

The  seemingly  complicated  formulation  of  the  second  postulate 
comes  from  the  fact  that  S is  only  defined  for  equilibrium  states, 
which  in  turn  were  defined  for  simple  systems.  For  our  purposes,  each 
species  (or  subspecies)  can  be  viewed  as  a simple  system  which  is 
initially  separated  by  walls  from  all  other  species.  The  composite 
system  is  then  in  an  internally  constrained  equilibrium.  When  all 
internal  constraints  are  removed  the  system  will  seek  a new  equilibrium 
state,  the  state  of  interest  in  this  work. 

3*6  Fundamental  Equations 

The  second  postulate  asserts  a functional  relation  between  S 
and  the  extensive  variables  U,  V,  x.  If  this  relation  is  known,  one 
can,  in  principle,  find  all  equilibrium  properties  of  the  system. 
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Such  a functional  relation  is  called  a fundamental  equation  (or 
characteristic  equation) . Since  the  free  energy  is  dervied  from 
fundamental  equations,  we  shall  now  study  some  of  their  properties. 


I 


Proposition  3-3:  Let  S = ¥(U,V,x).  Then  ¥ is  a homogeneous 

function  of  degree  one. 

Proof:  By  definition,  a function  f : V -» W is  homogeneous  of 

degree  p if  f(Ay)  = A^f(y)  for  all  A > 0 and  y € V.  Consider  N 


identical  simple  systems  in  equilibrium.  System  i is  characterized 
by  S*,  U*,  V1,  and  x*.  For  the  composite  system  formed  by  these  N 
simple  systems  we  have  by  postulate  2 (iii) 


S = ¥(U,V,x)  = ^(NU1,NVi,Nxi)  = N^U^V^x1) 


'v 


Extending  the  above  to  any  real  N we  obtain  the  result.  □ 

Proposition  5.4:  Let  S = ¥(U,V,x).  Then  ¥ is  invertible  with 

respect  to  U. 

Proof:  Clear  from  postulate  2,  (i)  (ii).  □ 

The  above  proposition  offers  another  way  to  write  fundamental 
equations,  with  U substituting  for  S as  the  dependent  variable. 
Thus  we  can  write  U = U(S,V,x). 


I 

; 

! 
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It  is  easily  verified  that  U is  also  a homogeneous  function 
of  degree  one.  This  property  of  U (and  as  we  shall  see,  also  of  the 


free  energy  F)  plays  a major  role  in  most  of  the  subsequent  parts 
of  this  work.  To  simplify  the  terminology,  we  shall  occasionally 
refer  to  "homogeneous  functions  of  degree  one"  simply  as  1-homogenous 
functions . 

Proposition  3.5: 


U = 


dU! 

dU 

dSl 

1 

S + 

V,x 

5v 

V t £ ||i-| 

S,x  J-1  \0xj/s,v,x1  J 


Proof:  By  the  previous  remarks  U = U(S,V,x)  is  1-homogeneous 

and  differentiable  (over  some  domain)  since  S is.  The  result  then 
follows  directly  from  the  well  known  Euler's  Theorem  which  states 
that  if  f(x)  is  differentiable  and  homogeneous  of  degree  p then 

p f ( x ) = x'Vf(x)  □ 


This  result  and  some  of  its  implications  were  discussed  by  Dantzig 
[15,17]. 


The  propositions  presented  here  set  the  stage  for  a similar 
treatment  of  free  energy  functions.  To  keep  the  discussion  in  the 
most  general  terms,  we  shall  derive  free  energy  functions  from  funda- 
mental equations  by  Legendre  transforms . 
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3.7  Legendre  Transforms  and  Homogeneous  Functions 


Let  G = G(x)  be  a differentiable  function  on  an  open  set 
Q c ln.  Let  y = y(x)  = \7G(x)  for  x £ G.  The  function 


cp(y)  = G(x)  - x'y  x £ a,  y = VG(x)  (l) 


is  called  the  Legendre  Transform  of  G(x). 

Let  H(x)  be  the  matrix  of  second  partial  derivatives  (assuming 
that  they  exist) 

dy. (x) 


By  the  implicit  function  theorem,  if  H is  nonsingular  there  exist 
functions  (2^,  i = 1,  2,  . . . , n such  that 

\ = ^(y) 

In  that  case  x can  be  eliminated  from  (l).  By  the  same  token  we 
could  write 

cp(x)  = G(x)  x'VG(x)  (2) 

or 

cp(x,y)  = G(x)  - x'y  (3) 


Legendre  transforms  (l)  replace  the  independent  variables  x 
by  the  partial  derivatives  of  the  function  G.  The  specific  form 
of  <p  will  be  (l),  (2)  or  (3)  depending  on  whether  y,  x or 
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both,  respectively,  serve  as  independent  variables.  In  the  latter 
case  only  part  of  the  x's  and  the  y's  are  independent  since  we 
assume  that  G has  at  most  n independent  variables. 


Lemma  3.6:  Let  f(x)  be  differentiable  on  an  open  set  G c ln. 

If  f(x)  is  homogeneous  of  degree  one  then  its  Legendre  transform 
is  cp(x)  = o. 


Proof : f(x)  = xVf(x)  by  Euler's  Theorem.  The  result 

follows . □ 


Lemma  3.7:  Let  f(x)  be  a homogeneous  function  (of  degree  one), 

differentiable  on  flc  En.  Then  y^x)  s dfCx)/!^  is  homogeneous 
of  degree  zero  for  i = 1,  2,  ...  , n. 


Proof » 


n 


f(Ax)  = Af(x)  = A £ y. (x)*x.  for  all  A > 0 

i=l  1 1 


but 


n n 

f(Ax)  = E y.(Ax)-Ax  = A E y.(Ax)-x 
i=l  i=l  1 


We  conclude 

n 


E y (Ax)-x 
i=l  1 1 


E y,  (x)-x 
i=l  1 1 


for  all  x € n,  A > 0 


Thus 


y^Ax)  = yi(x)  = iPy^x)  . 
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The  last  Lemma  shows  that  Legendre  transforms  can  be  used 
to  replace  extensive  variables  like  S,  V (which  are  homogeneous  of 
degree  one  since  they  are  proportional  to  the  mass  of  the  system) 
by  intensive  ones  which  are  homogeneous  of  degree  zero  and  thus  inde- 
pendent of  the  mass  of  the  system. 

We  can  thus  define  for  U = U(S,V,x) 


T = I rr 


dU 

dS 


(4) 


V,x 


P 3 - 


Idyl 

|PV, 


(.5) 


S,x 


“l  * 


dU 

^i/s,V,x 


i = 1,  2, 


(6) 


The  intensive  variables  T and  P defined  above  are  the  well  known 


absolute  temperature  and  pressure,  while  is  the  chemical  potential 


of  species  (see  section  2.4). 

The  complete  Legendre  transform  of  U is,  by  Lemma  J.6, 


cp(T,P,p)  = U - TS  + FV  - £ n.x  = 0 

i=l  1 


(7) 


The  quantity  U - TS  + FV  is  a partial  transform  with  respect 
to  S and  V.  We  shall  now  formalize  this  useful  notion  [11,  p.  96]. 


mil 
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Let  G(x)  be  a differentiable  function  on  ft  c if1.  Let 
O be  an  integer  subset  of  {1,  2,  ...  , n} . The  partial  Legendre 
transform  with  respect  to  the  variables  x^,  i £ a is  defined  as 


cpjx,y)  = G(x)  - E X y 
i £ a 


where  y £ A(g),  x £ ft.  Assume  without  loss  of  generality  that 
a = {1,  2,  ...  , k}  and  p = (k+1,  k+2,  ...  , n}  for  some  k. 
In  most  applications  of  (8)  we  choose  as  independent  variables 


sets  y^i  £ a)  and  x^.  (j  £ p) . We  write 


<p0u,y)  Vs’  •••  ’ Xn’  yi>  y2’  •••  ' 


Free  energy  functions  can  now  be  defined  in  terms  of 
partial  Legendre  transforms  of  U = U(S,V,x).  Let  a = (1,2} 
corresponding  to  S and  V.  Define 


F(x,T,P)  3 cpQ(x,T,P)  = U - TS  + FV  (9) 


It  follows  from  (7)  that 


F(x,T,P)  = E 

i=l  1 


^ = ^(xjTjP)  . 


mm 
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Equation  (lO)  is  identical  to  (2.4.1)  which  served  as  our 
original  definition  of  the  Gibbs  free  energy  function.  Similar  to 
F we  can  define  other  free  energy  functions  by  suitable  choice  of  a. 
Letting  a = ( 1)  we  obtain  the  Helmholz  free  energy  A, 

A(V,x,T)  a cpa(V,x,T)  = U - TS  (ll) 


and  by  (7),  the  Helmholz  free  energy  A becomet 


A(V, x, T)  = Z x.g.(V,x,T)  - V-P(V,x,T) 
i=l  1 1 


The  enthalpy  H is  defined  by  letting  a = [2) 


H(S,x,P)  = U + FV  = T(S,x,P)-S  + £ x g.(S,x,P) 

i=l 


(12) 


Our  main  goal  is  to  derive  general  properties  of  free  energy 
functions,  based  only  on  the  postulates  and  Legendre  transforms. 

Two  basic  properties  are  of  interest  in  our  discussion:  homogeneity 

and  convexity. 


Theorem  3-8:  Any  partial  Legendre  transform  of  a fundamental  equation 

is  a homogeneous  function  of  degree  one  in  the  extensive  variables. 


Proof i We  have  already  shown  that  a fundamental  equation  is 
homogeneous  of  degree  one  in  Proposition  3.3,  so  let  the  equation  be 
in  general  Xq  = G(X^,  X^,  •••  , X^).  By  Lemma  3.6  we  know  that 
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cp(X,Y)  = G(X)  - X'Y  = 0 


where  X = (X^  Xg,  . . . , Xn)  £ n c En,  Y = VG(x) . Let  a c (1,2, . . . ,n), 


P = (1,  2,  ...  , n}  - a.  Then 


cp  (X)  = G(X)  - E X Y (X)  = E X Y (X) 
i£a  J J 


where 


yj(x)  - ^ 


Therefore,  for  any  A > 0 


cPa(NO 


i 'V 


A E X.  = A E x Y (X) 

j£p  J A dxj  jep  J J 


%(X) 


This  result  is  not  entirely  new.  It  was  intuitively  known  and  proven 
in  the  past  (see  for  example  [17] ) for  specific  functions.  Our 
approach  is  a generalization  which  applies  to  all  thermodynamics  free 
energy  functions,  without  reference  to  their  specific  functional  form. 

The  foregoing  discussion  illuminates  the  reason  for  using 
F for  reactions  under  constant  pressure  and  temperature.  Under  these 
conditions  p,  (x,P,T)  = (^(x)  sc  that  F is  a function  of  composition 
only.  Under  constant  V and  T we  use  the  function  A (ll),  whence 
^(V,x,T)  = ^(x)  and  P(V,x,T)  = P(x),  and  again,  A is  a function 
of  composition  only. 
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The  homogeneity  of  F explains  the  reference  to  p^(x) 
as  the  "partial  molar  free  energy."  At  fixed  T and  P we  have 


F(x)  = £ M1(»)'X1  ■ l 


Moreover,  by  Lemma  3.7,  p^(x)  is  a homogeneous  function  of  degree 
zero  for  i = 1,  2,  ...  , n,  which  accounts  for  the  fact  that  in  most 
cases  ji(x)  = n(x). 

Many  difficulties  in  analysis  and  computation  of  free  energy 
functions  arise  from  their  bad  behavior  for  compositions  which  are 
not  strictly  positive,  i.e.,  which  have  vanishing  species.  We  observe 
that  right  continuity  of  F always  requires  that  F(0)  =0  since 


lim  F(Ax)  = lim  AF(x)  = 0 
1 — > 0 A -»  0 


Loomis  and  Steinberg  [37]  show  that  a 1-homogeneous  function  F is 
differentiable  at  0 if,  and  only  if,  it  is  an  affine  function.  Since 
the  free  energy  is  not  linear,  its  boundary  pathologies  are  in  this 
sense  a result  of  its  1-homogeneity. 


3.8  Convexity  of  Free  Energy  Functions 

The  postulates  of  Section  3-5  were  defined  in  terms  of 
the  entropy  S.  Tisza  £ 50]  and  Callen  [ll]  show  that  these  can 
be  described  equally  well  in  terms  of  the  internal  energy  U,  by 


I 


exchanging  the  roles  of  U and  S in  postulate  2 and  by  replacing 
the  "maximize  S"  in  (iv)  with  "minimize  U. " Internally  unconstrained 
equilibrium  is  then  reached  when  the  minimum  of  the  internal  energy 
is  attained. 

Based  on  postulate  2 we  can  now  prove 


Lemma  5,9:  U = U(S,V,x)  is  a convex  function. 


Proof;  Let  the  subscripts  c and  u denote  internally  con- 
strained and  unconstrained  systems  respectively. 

Consider  two  simple  systems  characterized  by  (S'jV'jX1)  and 
(s",V",x"),  respectively,  which  form  a composite  system.  From  postulate 
2 (iii)  we  have 


UG ( S ' , V ' , x ' ) + Uc(s",v",x")  = Uc(S  ' + s",  V'  + v",  x'  + x") 


By  the  energy  version  of  (iv) 


U (S'  + S",  V'  + V",  x'  + x")  > U (S'  + S",  V'  + V",  x'  + x") 
c u 


Since  we  are  dealing  only  with  equilibrium  states  we  can  now  eliminate 
the  subscripts  and  combine  to  get 


U(S',  V',  x* ) + U(S ",  V",  x")  > U(S 1 + s",  V'  + V",  x'  + x") 
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Homogeneity  of  U requires 


U(XS,  XV,  Xx)  = XU(S,V,x) 

Therefore 

U(XS\  XV',  Xx')  + U((1-X)S”,  (l-X)v’',  (i-X)x") 

= Xll(S  1 , V 1 , x 1 ) + ( 1-X)  U(s",V",x") 

> U(XS'  + (1-X)S",  XV'  + (l-X)V",  Xx'  + (l-X)x") 

This  completes  the  proof.  □ 

Theorem  3.10;  Let  fl  c fin,  and  let  G :fl  -»  S be  a twice  continuously 
differentiable,  homogeneous  function  of  degree  one  on  (J.  Let 

(x)  = G(x)  - £ y (x)-x 

i£a  1 

be  a partial  Legendre  transform  of  G with  respect  to  x^,  i £ a,  where 
h<x)  = i € N 3 (1 , 2,  ...  , n) 

Under  these  conditions,  cp  (x)  is  convex  in  the  untransformed" 

variables  x,,  j OL  for  all  a,  if  and  only  if,  G(x)  is  convex. 

J 
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Proof : The  "only  if"  is  trivially  proved  by  taking  a = 0 (the 


empty  set).  Then  cp^(x)  = G(x)  which  is  convex. 
Now  suppose  G is  convex,  and  let 


dy, (x) 

Vx)  • -ET 


SxTdxT  ’ 

i J 


i,  J € N 


The  symmetric  matrix  Q (x)  = (Q.  .(x))  must  be  positive  semidefinite 

^ J 

(see  Zangwill  [55>  P*  30])*  Accordingly,  all  principal  submatrices 
of  Q are  also  positive  semidefinite.  Let  a ^ N,  a / 0.  By  Lemma  3.6 


”a(x)  * 2 h(x)'xi 

i€N-a 

From  Theorem  3-8,  c^(x)  homogeneous  of  degree  one.  Hence 
^rx)  ‘ • 1 6 


Therefore 


d2y  (x) 

a 

dx.  dx . 
i J 


dy.  (x .) 
1 


dx.  dx . 
1 J 


Q_(x)  for  i,  j £ N-a 


The  second  partial  deriva-ives  of  9 (x)  with  respect  to  the  untrans- 
formed variables  form  a principal  submatrix  of  Q which  is  positive 
semidefinite.  Therefore  cp^(x)  is  convex.  □ 
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Theorem  3. 10  proves  that  free  energy  functions,  defined  as 
partial  transforms  of  fundamental  equations,  are  convex  in  the  extensive 
variables.  In  our  <‘asp 

F = F(T,P,x)  is  convex  in  x. 

Homogeneity,  on  the  other  hand,  shows  that  F cannot  be  strictly 
convex,  for  Let.  x"  = yx 1 for  some  positive  scalar  y.  Then 

F(Xx'  + (l-X)x")  --  F([X  + (l-X)y]x')  = [X  + (l-X)y]  F(x') 

= XF(x')  + (l-X)  F(yx')  - XF(x')  + (1-X)  F(x") 

Thus  F(x)  is  linear  for  x = x'  + yx' . 

We  have  shown  so  far  that  convexity  and  homogeneity  of  U 
and  F are  outcomes  of  the  basic  postulates.  We  shall  now  prove  the 
converse,  namely,  that  any  convex  homogeneous  function  will  satisfy 
the  minimum  energy  principle. 

Propositi'  Let  0(x)  be  a convex  homogeneous  function  defined  on 

n , ,. 

a convex  set  1 c 8 . Let  X and  X denote  respectively  the 

*,  , 

extensive  variables  of  two  simple  systems,  and  .let  G (X)  denote  the 
equilibrium  value  of  G(X).  Then 

G*(X’  + X")  _ G*(X 1 ) + G*(X")  . 
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Proof:  By  homogeneity  we  have  for  A > 0 


G*(X 1 + XX")  = (X+l)  G*  x'  + X" 


Together  with  convexity,  this  leads  to 


G*(X ' + XX")  < (X+l)  ^ G*(X')  + ~ G*(x") J 


= G*(X* ) + XG*(x") 


Letting  A = 1 we  get  the  result.  _ 

We  conclude  this  chapter  hy  noting  that  a Legendre  transform 
of  a function  can  be  geometrically  interpreted  as  the  envelope  of 
tangent  hyperplanes  to  the  graph  of  the  function.  As  such,  Legendre 
transforms  are  closely  related  to  "conjugate  functions"  which  play 
a major  role  in  Rockafellar 's  convex  analysis  [43].  When  a function 
is  replaced  by  its  partial  derivatives,  a constant  of  integration  must 
always  be  added  when  returning  to  the  original  function  since 
differentiation  "loses"  some  information.  This  is  the  reason  for 
cp(x)  being  defined  the  way  it  is  and  not  by  x'y*  The  last  point 
is  demonstrated  here  by  defining  the  chemical  equilibrium  problem  in 
terms  of  the  intensive  variables  for  an  ideal  system  under  fixed 

temperature  and  pressure. 
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We  have 


Hj  = ^j(x)  = Cj  + log(x  /^J  , 

Hence 

*j  * V^J  ' “V  ’ 

n 

- £ MA(j)  ' cj> 


j £ 00 


j € <k>  (1) 

(2) 


The  mass  balance  constraints  are 

Z x,  Z a exp(|i  -c  ) = b,  i=l,2,  ...  , m (3) 

k ^ J£(k)  1J  J J 1 

The  variables  x^  are  in  essence  the  "integration"  constants  for 
the  free  energy  of  each  phase. 
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GEOMETRIC  PROGRAMMING  AND  CHEMICAL  DUALITY 

The  first  part  of  this  chapter  reviews  the  duality  theory  of 
geometric  programming  due  to  Duffin,  Peterson  and  Zener  [25] • The 
review,  including  some  important  theorems,  is  followed  by  a dual 
formulation  of  the  chemical  equilibrium  problem.  We  shall  look  at 
several  forms  of  the  dual  chemical  problem  and  then  try  to  give  a 
chemical  interpretation  relating  these  problems  to  well  known  chemical 
concepts  and  laws.  In  this  chapter  we  shall  at  times  sacrifice 
mathematical  rigor  for  brevity. 


4.1  Primal  Geometric  Programs 

In  the  sections  dealing  with  geometric  programming  we  shall 
rely  mostly  on  Duffin,  Peterson  and  Zener's  notation  [25],  with  some 
changes  to  illuminate  the  relation  to  chemical  problems.  The  theory 
of  geometric  programming  was  developed  primarily  for  its  engineering 
applic  ■‘‘ions,  which  take  the  form  of  the  following  primal  geometric 
program  t 


kk 


Program  PGP;  Find  a vector  t € fi  so  as  to 

Minimize  gdt.)  = Z u.(t) 

0 je<o)  J 


•d) 


Subject  to  g,  (t/ 


Z u>t)<i,  k - 1,  2 , 
j€<k>  -J 


and 


t H 


(V  V “ ' Xj 


> 0 


. , K (2) 


(3) 


Here 


m a 


u,(t) 

J 


r n t iJ 

u i i ’ 

d i=l 


j = 1,  2, 


C.  are  positive  real  numbers,  a are  real  numbers,  the  sets  (k) 
J 1 J 

for  k = 0,  1 y 2y  ...  , K are  integer  sets  partitioning  the  set 
N = (1,  2,  . . . , n}  as  described  in  Section  2.2.  For  convenience 
we  also  define  the  integer  sets 


M = { 1,  2,  ...  , m} 
K = U,  2,  ...  , K} 


The  terms  u.(t)  resemble  positive  general  polynomial  terms  and  are 
J 

therefore  called  posynomial  te  Similarly,  the  functions  g.  ft), 

k £ { 0}  U K are  called  posynomial  functions  or  simply  posynomia  ls . 

The  matrix  A = {a.  .)  is  the  exponents  matrix,  and  the  positive 
vector  C = (Cj_,  Cg,  ...  , Cn)  is  the  coefficients  vector.  Problem 


PGP  can  therefore  be  completely  characterizted  by 


UfWH-  IM'WJ  I.  JJWIWUiS.IJI'  1 


(i)  An  m X n exponento  matrix  A 

(ii)  A positive  n-vector  of  coefficients  C 

(iii)  A partition  of  the  set  N defining  the  sets  (k),  k = 0,  1,  K. 

We  say  that  program  PGP  is  consistent  if  there  exists  a vector  t satisfy- 
ing (2)  and  (3).  The  program  is  super-consistent  if  there  is  a t > 0 
for  which  (2  ) are  all  satisfied  with  strict  inequality. 


4.2  Dual  Geometric  Programs 

With  each  primal  geometric  program  there  is  an  associated 
dual  geometric  program. 


Program  PGP;  Find  vectors  8 € En  and  A 6 EK  so  as  to 


T n B,-j  r K \ -1 

Maximize  v(&)  = II  ( c j ) [ n \ J 

(1) 

Subject  to  Z 8 - -1  (Normality  condition) 

je<o>  3 

(2) 

A8=  0 (Orthogonality  conditions) 

(3) 

6 = (5^,  Sgj  •••  ) 6 ) > 0 (Nonnegativity) 

(4) 

and 

\ = E B.,  k e k 

* J€(k)  J 

(5) 
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Here  the  coefficients  C , the  matrix  A and  the  partitions  (k> 

J 

of  N are  the  came  as  defined  for  FGP.  Tc  maintain  continuity  of 
v(&)  on  the  nonnegative  orthant.  of  R'1  we  define 

c 

5 - b 2 0 when  5=0 

Note  that  the  dual  program  an  he  completely  characterized  by  the  same 
information  characterizing  PGP.  Frogram  DGP  is  said  tc  be  consistent 
if  there  exis-ts  a 5 > 0 satisfying  (2)  and  .(,5). 


U . 5 Duality  Tneory  of  Geometric  Programming 

Tht  mail  1 emma  of  geometric  programming  states  the  relation 
between  the  primal  and  dual  programs. 


Lemma  4.1:  (Duff in  et  al.  [2j,  p.  llU] ).  If  t satisfies  the  con- 

straints of  priinaL  program  PGP,  and  5 satisfies  the  constraints  of 
its  dun'  program  DGP  then 

g (t)  > v(5)  (l) 


Moreover,  under  the  same  conditions  g^(t)  = v(b)  if,  and  only  if, 


u,(t)/g  (t)  j 6 (0) 

d o 


(?) 
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J 


j e (k),  k € K 


The  lemma  serves  as  a basis  for  the  following  duality  theorem: 

Theorem  4.2:  ([25,  p.  80]).  Suppose  primal  program  PGP  is  superconsi stent, 

and  gQ(t)  attains  a minimum  value  at  a feasible  point  t*.  Then 

(i)  Dual  program  DGP  is  consistent  and  attains  its  maximum  at  a 

* 

dual  feasible  point  6 . 

(ii)  v(6*)  = gQ(t*) 

(iii)  There  exist  nonnegative  Lagrange  multipliers  tj,  , k € K such 

K, 

that  the  Lagrange  function 

K 

L(t,T))  E gQ(t)  + E \tgjt)  - 1]  (5) 

k=l 

has  the  property 

L(t  >*))  < 60(t  ) = L(t  ,T)  ) < L(t,T)  ) (4) 

for  arbitrary  t > 0 and  arbitrary  r\  > 0.  Moreover)  there 

* 

exists  a maximizing  vector  6 for  DGP  such  that 

UjU'  > J € (o) 

(5) 

\uj(t  )/60(t*)  » j € (k>  , k € ic 

Furthermore, 


V8  > - \/60(t  > 


k € K 


(6) 


(iv)  If  5 is  a maximizing  point  for  dual-  program  DGP,  each 

* 

minimizing  point  t for  PGP  satisfies  the  system  of  equations 


»,(*•>  * 


6^  v(b*) 


j £ (o) 

j (k) 


(7) 


where  k 


ranges  over  the  integers  for  which 


\(M  > 0. 


A dual  program  is  said  to  be  canonical  if  there  exists  a 
positive  vector  b > 0 satisfying  the  dual  constraints  (2.2)  and 
(2.5).  Otherwise  the  program  is  said  to  be  degenerate . We  shall  say 
that  a primal  program  is  canonical  or  degenerate  when  the  corresponding 
dual  program  is. 


Theorem  h . ([25,  p.  169] ).  Suppose  dual  program  PGP  and  its  corre- 
sponding PGP  are  canonical.  Then  DGP  is  always  consistent,  but  PGP  is 
consistent  if,  and  only  if,  DGP  has  a finite  positive  maximum.  More- 
over, under  these  conditions  the  minimum  of  RIP  is  attained  at  some 

if 

finite  primal  feasible  point  t >0  and  is  equal  to  the  maximum  of 
DGP. 
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4.4  Transformed  Primal  Programs 


For  theoretical  and  practical  reasons,  it  is  sometimes 
advantageous  to  consider  a transformed  version  of  PGP.  Let 


zi  = log  ti 


i = 1,  2,  ...  , m 


The  transformed  (primal)  geometric  program  is 


Program  TGP 


Minimize  gn(z)  - E u (z) 

U j£  (o)  J 


where 


Subject  to  g.  (z)  2 E u (z)  <1,  k€  R 

j£<k>  J 


u.(z)  = exp(-c  + E a z } , j £ N 


=J  - - log  Cj 


We  used  identical  notation  for  PGP.  It  will  always  be  clear  from 
context  whether  we  are  considering  program  PGP  or  TGP. 

Two  points  should  be  noted  about  program  TGP 
(i)  The  variables  z i are  unrestricted  in  sign. 

(ii)  The  program  is  convex  since  u.(z)  are  convex  functions. 


4.5  Chemical  Duality 


Following  the  development  of  the  optimality  conditions  by  Kuhn 
and  Tucker  [56],  it  was  shown  by  White  et  al.  [54]  and  later  by  Dorn 
C25l  that  the  Lagrange  multipliers  of  the  chemical  equilibrium  problem 
CPi  give  rise  to  a new  problem  called  a dual  chemical  problem.  To 
avoid  confusion  between  geometric  programming  duali-y  and  chemical 
duality  we  shall  refer  to  problems  CPI  and  CPN,  i.e.,  to  the  chemical 
equilibrium  problems,  as  primal  chemical  problems.  Dual  chemical 
problems  are  those  problems  involving  the  multipliers  of  problems  CPI 
or  CPN  of  Section  2.5. 

White's  and  Dorn's  dual  chemical  programs  are  not  "pure"  duals 
in  that  they  contain,  in  addition  to  multipliers,  also  primal  chemical 
variables,  namely,  composition  variables.  Avriel  [4],  and  Pussy  and 
Wilde  [40],  showed  that  geometric  programming  can  be  applied  to  generate 
a pure  dual  chemical  problem,  in  fact, they  showed  that  the  chemical 
equilibrium  problem  is  equivalent  to  the  dual  of  a geome^ri  program 
(DGP).  It  is  then  a straightforward  task  to  show  that  a dual  chemical 
problem  equivalent  to  PGP  can  be  formulated. 

To  show  that  problem  CPJ  is  equivalent  to  DGP  we  define 

— ~bj  i = 1,  2,  ...  , m 

bQ  = = 1 and  <0)  = ID) 

We  then  identify 
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so  that 

Xj  " 8jXk(j) 

Rather  than  minimize  F(x)  in  (2.5.1),  we  can  maximize  exp(-F(x)), 
which  after  the  above  transformations  is  seen  to  be  v(b).  The 
orthogonality  conditions  A&  = 0 are  satisfied  for  the  augmented  matrix 
A which  has  the  added  column  Aq  = -b. 

There  is  no  difficulty  in  showing  the  converse--that  DGP  can  be 
formulated  as  a CPI  so  that  the  relation  is  indeed  two-sided.  We  remark 
here  that  the  general  DGP  is  equivalent  to  an  abstract  chemical  problem. 
The  special  properties  of  matrix  A in  chemical  problems  (for  example, 
Proposition  3.1)  do  not  hold  for  a general  DGP. 

Our  interest  will  be  focused  mostly  on  the  dual  chemical  problem. 
In  light  of  ihe  transformations  above,  it  has  the  following  form: 

Problem  DCP  (Dual  Chemical  Problem) 

Minimize  -log  g^(t)  H *b ' log  t £ - £ b^  log 

Subject  to  g,  (t)  = Z u (t)  <1  k £ K 

k J£<k>  J " 


t,  > 0,  t_  > 0,  ...  , t > 0 


Here 


m a.  . 
u.(t)  = C.  II  t.1J 
J JH  1 


j c.  N 


» exp(-c^) 


j vJ  N 


where  c.  are  the  free  energy  coefficients f A = {a. .}  and 
J i J 

b = (b,,  b_,  , b ) are  the  matrix  and  right  hand  side  of  the 

'12  m 

balance  equations. 


mass 


log  t = (log  t±,  log  tg,  ...  , log  tm) 


Another  formulation  of  the  chemical  dual  is  based  on  the 
transformed  geometric  program.  This  form  is  essentially  the  one  used 
by  Dorn  [2jj]  and  Bigelow  7],  and  it  simply  replaces  log  t^  by  z.^ 
so  that  we 


where 


Minimize  -leg  g q(z)  = *b'z 

Subject  to  g (z)  = Z u.(z)  <_  1 

j€(k) J 


k t K 


m 


.(z)  = exp(-c .+  Z a.  .z.  ) j t i 


J U1  U i 
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k . 6 Interpreting  +hf.  Eual  Chemical  Problem 

Most  of  the  implications  and  applications  of  cnemical  duality, 
especially  those  of  equivalence  to  geometric  programs  will  be  treated 
in  the  following  chapters.  In  the  rest  of  this  chapter  we  shall  study 
the  physical  meaning  of  the  dual  chemical  problem.  Specifically,  we 
shall  interpret  the  dual  chemical  variables  and  constraints. 

Dorn  [23]  and  Whi*e  et  ad  who  realized  the  existence  of 

duality  in  chemical  equilibrium,  also  attempted  to  interpret  the  dual 
variables.  It  is  easily  seen  that  dual  chemical  variables  t or  z 
relate  to  the  subspecies  and  treir  mass  balance  equations.  Tne  form 
of  the  objective  function  in  DTP  (which  equals  F(x)  at  an  optimal 
t and  x ) suggests  that  log  t.  is  some  energy  measure  of  subspecies 
II  • Indeed  all  the  interpretations,  including  a later  one  by  Duffin 
and  Zener  [26]  which  considered  entropy  functions,  followed  this  line. 

We  find  two  major  defects  in  all  these  attempts.  First  none 
of  the  above  approaches  accounts  for  the  facT  that  the  dual  chemical 
variables  for  a giver,  chemical  equilibrium  problem  are  not  unique. 
Second,  all  are  quite  'unintuitive  and  do  not  readily  relate  to  known 
chemical  concepts. 

Dorn  suggests  that  "*he  dual  variables  are  dimensionless 

energies,  contributed  by  the  elements  (subspecies)."  White  et  al. 

state  that  "it.  (the  multiplier  of  [b . - £n  , a.  ,x  ] ) is  the  free 
1 1 J=1  ij  j 

energy  contribution  due  to  the  presence  of  one  mole  of  atom  (sub- 
species) i."  The  set  of  subspecies  for  a given  system  need  not  be 


5^ 


l 


unique  (see  Section  2,j5).  Moreover,  the  free  energy  function  itself 
is  only  a relative  measure.  Its  value  depends  on  the  definition  of 
certain  standard  states  and  energies.  We  shall  deal  first  with  these 
states  and  then  return  to  interpret  the  dual  chemical  variables. 


4.7  Equivalent  Formulations  of  the  Chemical  System 
Consider  the  problem 

Problem  CPI 

Minimize  F(x) 

Subject  to  Ax  = b,  x > 0 

Where  x € £n,  A € if'1*0,  b £ fT . 

Let  D € be  a real  nonsingular  matrix  and  let  y 6 R™  be  a given 

vector.  Let  A = DA,  b = Db,  and  consider  the  following  problems: 

prob 1 em  CPC 

Minimize  F(x) 

Subject  to  Ax  = b,  x > 0 

Problem  CP3 

Minimize  F(x)  = F(x)  - y'Ax 
Subject  to  Ax  = b,  x > 0 
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It  is  easy  to  see  that  problems  CPI,  CP2,  and  CP3  are  equivalent 
in  the  sense  that  x is  a solution  to  one  of  them  if,  and  only  if,  it 
is  a solution  to  all  of  them.  The  equivalence  of  CPI  and  CP3  was  used 
by  Dantzig  et  al.  [l8]  to  obtain  more  convenient  coefficients  in  F(x). 
In  the  ideal  case 


F'(x)  = x'  (c  + log  x) 


Hence,  the  objective  function  of  CP3  becomes 


F(x)  = x'-(c  - A'y  + log  x) 


Looking  at  the  dual  chemical  problems  associated  with  CPI, 

CP2  and  CP3  we  see  that  for  CPi  (in  the  ideal  case)  the  dual  chemical 
problem  is  DCP. 

For  CP2  we  have  a similar  dual,  with  a.  . replacing  a . and 

^ J lj 

b^  replacing  b^ . For  CP3  we  obtain  the  dual  chemical  problem  by  re- 
placing C,  of  DCP  witn 

J 


C . = exp[-c  , 4y'A  . 1 
J J J 


* * 

Obviously  the  solutions  t or  z to  these  dual  chemical  problems 
are  not  the  same. 

In  a more  general  context,  let 


F(x)  - Z |i  ,(x)  -x  . 
j=l  J J 

which  is  the  most  general  form,  as  shown  in  (2.4.2),  (3.7.I3), 
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Kuhn  Tucker  Conditions  applied  to  these  problems  then  state 

that  if  F is  differentiable  and  a constraint  qualification  holds, 

# 

a necessary  condition  for  x to  be  a solution  to  CPI,  CP2  or  CP3, 
respectively,  is  that  there  exist  vectors  z € such  that 


For  CPI 


z*A  = u'(x*) 


For  CP2 


z A = p (x  ) 


For  CP3 


z A = p'(x  ) - y'A 


Since  the  values  of  z in  (5)>  (6)  and  (?)  are  clearly  not 
the  same,  the  question  arises  as  to  what  is  the  nature  of  the  differ- 
ences between  the  original  chemical  problem  CPI  and  the  two  equivalent 
problems  CP2  and  CP3- 

Problem  CP2  amounts  to  a redefinition  of  the  subspecies.  The 
only  restrictions  on  the  subspecies  being  that  they  be  independent  and 
sufficient  to  describe  all  the  species  uniquely,  one  can  apply  non- 
singular linear  transformations  on  the  subspecies  and  maintain  the 
conditions  above  (although  the  resulting  set  of  new  subspecies  may 
appear  strange  to  a chemist).  We  demonstrate  such  a transformation 
by  the  following  example. 


Example  4.4.  Consider  a single  liquid  phase  system  with  the  set  S 
of  species  containing  COg,  H+,  OH  , HgO,  Ac,  HAc,  HgCOy  C0“  where 
Ac  represents  acetic  ion.  The  set  B of  subspecies  is  composed 
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Here  the  new  column  of  subspecies  is  obtained  from  the  previous  one 
by  premultiplying  it  by  (D 


1 0| 


0 0 


C02  - H + OH 


1 0 


0 1 


H+  + Ac 


Thus,  linear  transformations  of  A amount  to  changing  the  set  of  sub- 
species B.  Notice  that  the  free  energy  coefficients  were  not  changed. 
How  can  one  change  the  subspecies  and  maintain  the  same  free  energy 
function?  The  answer  is  found  in  the  definition  of  the  free  energy 

coefficients  c . The  same  question  arises  with  regard  to  CP}.  In 
J 

this  case  we  can  write 


where 


F(x)  = x(c(y ) + log  x) 


c(y)  = c - A'y 


For  a given  fixed  vector  y,  problem  CP3  amounts  to  changing  only 
the  free  energy  coefficients. 

To  see  the  meaning  of  these  transformations,  we  review  the 

derivation  of  the  vector  c [19,  21].  To  compute  all  c with  a 

J 

common  reference,  it  is  assumed  that  each  species  S.  is  formed  from 

J 


the  subspecies  via  the  formation  reaction 


m 


Z a.  .B.  = S 


i=l 


ij  i 


(9) 


where  each  of  the  subspecies  is  in  its  reference  state.  A reference 
state  is  usually  chosen  to  be  the  state  where  the  pure  subspecies 
is  in  its  most  stable  form  at  25°C.  At  any  rate,  each  subspecies  has 
some  reference  state  and  some  associated  reference  free  energy  f . 
The  reaction  (9)  has  some  associated  change  in  free  energy,  since  it 
can  be  viewed  as  an  internally  unconstrained  equilibrium  (see  Section 
Let  this  change  be  y then 


m 


= + Z a.  .f 


'J  J 


i=l 


iji 


(10) 


where  c,  represents  the  dimensionless  free  energy  of  one  (pure) 

J 

mole  of  S..  (AF.  and  f.  are  assumed  dimensionless.)  Now  suppose 
J J 

that  the  reference  states  are  redefined  so  that  the  new  reference 
free  energy  of  B^  is  f^  + y_^ . The  new  coefficient  is 


m 


c . 
0 


+ £ “i/h  + 


= c . + 


m 

Z a 
i=l 


ij*i 


Thus 


c =•  c + A'y 


Problem  CP3  can  thus  be  interpreted  as  a system  with  different  reference 
states,  which  have  new  reference  free  energies.  Similarly,  the  fact  that 
c was  unchanged  in  CP2  is  explained  by  an  implicit  change  in  reference 
states,  associated  with  the  new  set  of  subspecies. 
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4.8  The  Physical  Significance  of  the  Dual  Chemical  Problem 


The  relation  of  the  dual  chemical  problem  to  geometric  pro- 
gramming and  our  observations  about  reference  states  will  now  be 
employed  to  interpret  the  variables  and  constraints  of  problem  DCP. 
Let  t*  be  a solution  to  DCP  (Section  4.5).  Let  the  primal  chemical 
problem  be  CPI  with  F(x)  given  by  (7-1)  (this  is  problem  CPI  of 

Chapter  2--the  standard  ideal  chemical  equilibrium  problem).  Let  x 

* 

be  a solution,  namely,  x is  an  equilibrium  composition  vector. 

Applying  the  equivalence  to  geometric  programming  and  the 
duality  theorem  4.2,  we  have  by  (4.3.7) 


u u jl  u ” c in  w o. , 

h = \ uj(t  > - 4 e J (V 


where 


4 = k(j) 


In  other  words,  when 


Xk(j)>° 


x . 
J 


m 

n 

1=1 


-c 


e 


j 


(1) 


This  relation  has  exactly  the  form  of  the  mass  action  law  (Section  3.2) 

for  the  formation  reaction  of  S.  (Equation  (7-9))>  where  K(0)  = K.  = e 

J J 

Indeed  we  know  from  classical  thermodynamics  (Denbigh  [21,  p.  140] ) 

that  the  (dimensional)  free  energy  of  S is  given  by  -RT  log  K . 

J J 

From  (7.10)  we  then  have 


6l 


RT  c = -RT  log  K. 
J J 


<2) 


which  then  reduces  ( 1)  to 


m 


x.  n (t.) 

J i=i  1 


*' "ai j 


(3) 


The  mass  action  law  (5)  stated  here  relates  concentrations 
(mole  fractions)  of  reactants  and  products  in  the  formation  reaction 
to  the  equilibrium  constant  of  the  reaction,  where  the  reactants 
(subspecies)  are  in  their  reference  states.  We  conclude  therefore 
that  the  dual  chemical  variables  t^  represent  concentrations  of  the 
subspecies  in  their  reference  states ■ 

This  interpretation  is  not  merely  an  abstract  construction. 

When  one  or  more  phases  of  the  system  do  actually  correspond  to  a 
reference  state  of  some  subspecies,  say  phase  is  the  reference 

K. 

* 

state  of  Bi,  then  the  cor  entration  t will  be  the  actual  equilibrium 
concentration  of  in  that  phase.  Note  that  we  can  always  assume 

without  loss  of  generality  that  is  included  also  as  a species 

in  each  phase.  In  practice  (depending  on  the  reference  states  chosen), 

* 

t^  is  usually  a very  small  number,  indicating  that  B^  does  not 

exist  free  in  the  system  in  any  significant  amount.  (See  Section  6.2.) 

Equation  (l)  yields  an  easy  interpretation  to  the  posynomial 

terms  u.(t)  of  the  dual  chemical  problem  DCP.  We  already  saw  their 
J 

relation  to  mass  action  laws  of  formation  reactions 
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Jk  'ft 


•'ililiHiWuittri 


e 


-c 


u.(t  ) = 
J 


m / *\ai  i ,# 
n t.  ^ x 

i=l  1 J 


(4) 


Hence,  at  equilibrium  each  posyriomial  term  represents  the  equilibrium 
concentration  (mole  fraction ) of  an  associated  species . A point  t 

which  is  feasible  but  not  optimal  represents  a state  in  which  the 

. > .* 

formation  reactions  are  not  completed  so  that  u.(t;  < x,.  Finally, 

J J 

the  dual  chemical  constraints 


g,  (t)  = £ u (t)  < 1 

* J€(k)  J 

are  a restatement  of  the  fact  that  mole  fractions  sum  to  unity  for 
each  existing  phase  at  equilibrium. 


Equations  (3.7)  imply  that  g (t 


1 whenever  xfc  > 0. 


Duff in  and  Zener  [26]  noted  that  for  a gas  phase  the  dual  chemical 
constraint  can  be  reformulated 


Z p.(t ) < p 

jc(k)  J 


This  is  Dalton’s  Law  which  states  tnat  the  partial  pressures 
of  ideal  gases  at  equilibrium  sum  to  the  total  pressure  P when 
the  gas  pnase  exists. 

The  basic  difference  between  the  primal  cnemical  variables 
x arid  the  dual  variables  t (or  z)  is  that  while  the  former  are 
extensive  variables  whose  values  depend  on  the  mass  of  the  system, 


1 


the  latter  are  intensive  variables  and  can  thus  be  viewed  as  the 


outcome  of  a Legendre  transform  of  an  extensive  function.  The  hidden 
"integration  constants"  required  to  recover  extensive  variables  from 


intensive  ones  are  the  multipliers  tj  (of  the  dual  chemical  con- 

K 

■ft  * 

straints),  which  are  needed  to  compute  6 {or  x ) in  equation  (3.5) 

J J 


6U 


CHAPTER  5 


PROPERTIES  OF  EQUILIBRIUM  SOLUTIONS 
5.1  Introduction  and  Notation 

In  this  chapter  we  shall  study  the  solutions  to  the  chemical 
equilibrium  problem.  The  main  issues  are  existence,  uniqueness, 
bounds,  common  properties  of  solution  points,  and  general  properties 
of  solution  sets,  henceforth  called  equilibrium  sets . These  points 
will  be  reviewed  in  light  of  previous  results  relating  to  Legendre 
transforms  and  geometric  programming  duality  theory. 

Credit  for  much  of  the  original  work  on  this  siibject  is  due 
to  Shapiro  and  Shapley  [47],  and  to  Bigelow  [7].  Some  of  our  results 
are  not  new,  and  were  elaborately  proved  in  [47].  We  shall  apply 
geometric  programming  duality  to  obtain  simpler  proofs  wnich  are 
direct  consequences  of  the  theory  developed  by  Duffin,  Pe+erson  and 
Zener  [25] • Other  results  here  generalize  previous  results  in 
two  ways:  first,  by  considering  more  general  free  energy  functions, 

and  second,  by  obtaining  results  under  somewhat  weaker  conditions. 

This  chapter  provides  a theoretical  basis  to  the  study  of  equilibrium 
sets  and  for  bridging  the  gap  between  the  mathematical  implications 
of  the  model  and  the  behavior  of  real  chemical  systems. 
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V, 


Unless  stated  otherwise,  we  assume  a closed  system  under 
constant  temperature  and  pressure.  The  system  is  defined  by  the 
triplet  (F,A,b}  as  described  in  Section  2.5,  where  F is  taken 
to  be 


F(x) 


n 

z 

j=l 


(1) 


When  F has  the  form 


F<x) 


n 

z 

j=l 


xJ(cJ 


log  x^) 


(2) 


we  refer  to  it  (and  to  the  system)  as  ideal.  F is  assumed  differ- 
entiable on  the  positive  orthant  of  E*1.  The  feasible  set  of  the 
system  is  denoted  by  X 

X = X(A,b)  s (x  € Bn|Ax  = b,  x > 0}  (5) 


Its  intersection  with  the  positive  orthant  E^  is 

X+  s lx  6 X jx  > 0}  . 

We  say  that  the  system  is  canonical  when  X+  / 0 where  0 denotes 
the  empty  set.  Otherwise  the  system  is  degenerate . (See  Section  4.3.) 
The  equilibrium  set  (solution  set)  of  a given  system  is  denoted 

D(F|X)  = (y  € X|F(y)  < F(x)  for  all  x £ X)  (4) 
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5.2  Some  Equilibrium  Set  Characteristics 

Theorem  5-1:  For  a real  chemical  system  (F,A,b),  the  solution  set 

C(F|X)  is  bounded. 

Proof;  By  Proposition  3-1,  X is  bounded  so  the  result  is 
immediate.  □ 


F was  shown  to  be  convex  on  its  domain  of  definition,  which 
is  normally  E^.  When  possible,  the  domain  is  extended  to  the  non- 
negative orthant.  The  following  simple  lemma  insures  that  in  these 
cases  convexity  is  maintained  on  the  extended  domain. 

Lemma  5.2:  Let  F : Q — > E be  a convex  function  on  the  open  set 


Let  si  denote  the  closure  of  En.  Suppose  that 

lim  F(x  + td)  exists  for  all  x and  0 in  En  for  which 
t ->  0 

x + th  £ Q for  all  0 < t ^ tQ 

Then  F is  convex  on  f: 

Proof;  Let  x,  y£  En.  Let  9 > 0,  so  that  x + te  £ 2 
for  all  t > 0.  Then,  for  any  X £ [0,1]  we  have 

FrA(x  + tO)  + (1-X) (y  + t0  ] < XF(x  + to)  + (l-A)  F(y  + to) 
Hence 


4 

4 2 
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F[Ax  + (l-A)y  + tt)]  < AF(x  + t0)  + (l-A)  F(y  + t0) 

Taking  the  limit  as  t > 0+  on  both  sides,  we  obtain  the  result.  □ 


Theorem  5-3:  If  F satisfies  the  conditions  of  the  preceding  lemma, 

then  D(Fjx)  is  convex. 

Proof:  D is  simply  the  minimum  set  of  a convex  program.  □ 

This  is  Lemma  9-3  in  [47],  applied  here  to  a general  free 
energy  function.  Since  convexity  on  K+  is  assured  by  Theorem  3-10, 
only  the  right  continuity  of  F on  X is  required  to  satisfy 
Lemma  5*2. 

Two  important  characteristics  of  solutions  which  were  developed 
via  "quasi  dependence"  in  [4-7]  are  shown  here  to  be  direct  results  of 
geometric  programming  duality. 


Tneorem  5.4:  Let  the  system  {F,A,b}  be  canonical,  where  F is  ideal. 

Let  x*  £ c(FjX).  Then 

* -* 

x . = 0 if  and  only  if  x,  . . » = 0 
J k(J) 


Proof:  Let  k = k(j). 

If  x^  = 0 clearly  x^  = 0 by  definition  of  x^.  Now 
suppose  Xj  = 0 but  x^  > 0. 

By  Theorem  4 >3  there  exists  a finite  positive  t solving 
the  dual  chemical  problem  DCP.  According  to  Theorem  4.2 
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Uj(t  ) = e J n (t^  iJ  = Xj/^  = 0 

# 

contrary  to  the  positivity  and  finiteness  of  t . □ 

This  well-known  result  has  far  reaching  chemical  implications. 
It  says  that  at  equilibrium,  either  a whole  phase  vanishes,  or  each 
of  the  species  in  the  phase  must  be  present  in  some  positive  amount, 
when  the  system  is  canonical.  This  justifies  our  remark  about 
equilibrium  concentrations  of  subspecies.  These  can  always  be  assumed 
to  exist  if  a reference  phase  exists  at  equilibrium  (Section  4.8). 

Theorem  5.5:  Suppose  (F,A,b)  is  canonical  and  F ideal. 

Let  x £ d(f|x)  and  y £ j)(F|x)  be  two  distinct  equilibrium 
solutions.  Then 


whenever  both  are  defined. 

Proof:  Theorem  4.3  guarantees  the  existence  of  finite  and 

positive  vectors  t(x)  and  t(y)  solving  the  dual  chemical 
problem  DCP.  Furthermore,  every  solution  t to  DCP  must  satisfy 
condition  (iv)  of  Theorem  4.2,  namely 

uj(t)'xA(j>  when  \(j)>0 

similarly 

uj(t)  = yA(j)  ',hen  yk(j )>0 
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Both  equations  must  hold  for  every  t,  therefore 


xA(j>  'yA(j> 


whenever  both  a~e  defined.  □ 


To  the  chemist,  this  result  comes  as  no  surprise.  It  indicates 
a certain  concentration  invariance,  i.e.,  some  uniqueness  of  concen- 
trations in  equilibrium  solutions.  It  may  be  harder,  though,  to  accept 
the  fact  that  the  same  phase  may  vanish  in  one  solution  but  not  in 
another  solution  to  the  same  system!  To  realize  such  a situation, 
observe  that  the  mathematical  formulation  does  not  prevent  us  from 
"splitting"  a phase  in  our  model  into  two  identical  phases  and  treat- 
ing these  as  two  completely  separate  phases. 

Suppose  is  a solution  for  (j  c k)  when  is  treated 

as  a single  phase,  that  is,  x.  (j  c k)  are  equilibrium  values.  Assume 

J 

x,  > 0 so  that  there  are  equilibrium  mole  fractions  x.,  j £ k. 

K J 

We  now  separate  into  ® and  ® and  construct 

k kx  kg 

a solution  to  the  new  system  by  letting  a £ [0,1]  and 


\ = 


\ ■ 


1 

x . = ax . 
J J 


x.  = [i-a]x. 

J J 
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Clearly  the  new  solution  satisfies  the  constrain*.  -5  ~4’  be  the 

/>1  /s2 

new  free  energy.  Since  x = x = x , we  obtain 

J J J 


F(x)  - F(x) 


Z x (c  + log  x ) - ( Z x^(c  + log  xb  + Z xbc  + log  x^)} 
je(k)  J J J€<k)  j J je(k>  J J J 

Z [x  - (x^  + xb]  (c  + log  X ) =0 
j€(k>  J J J J J 


This  implies  that  the  new  system  is  also  in  equilibrium, 
because  its  free  energy  cannot  decrease  below  F(x),  the  new  system 
being  "internally  constrained"  by  the  separation  of  Convexity 

of  o then  implies  that  for  all  A 6 [0,1] 


y = Ax^  + (l-A)x^  for  all  j 6 k 
J J J 


is  also  a solution,  since  0!  was  arbitrary.  In  particular,  A = 0 


and  A = 1 cause  4 and  4 respectively  to  vanish.  In  a 

k2  kx 

well  formulated  system  such  a situation  should  not  arise.  Conditions 


which  insure  that  phases  are  not  "split"  in  this  way  are  part  of 
the  regularity  condition  discussed  in  the  following  sections. 


5-3  Reduction  of  Degenerate  Systems 


A system  {F,A,b}  may  be  canonical  and  yet  the  solution  need 
not  be  strictly  positive  as  was  just  shown.  A solution  x is  degenerate 
if 

x € c(f|x;  but  x £ X+ 


This  section  deals  both  with  system  degeneracy  and  solution  degeneracy. 
In  both  cases  it  is  shown  that  the  problem  can  be  reduced  by  eliminating 
vanishing  species. 

Lemma  5.6:  A feasible  system  is  degnerate  if  and  only  if  there  exists 

an  index  set  Ic  N = 11,  2,  ...  , n),  I / 0,  such  that 
(i)  j £ I implies  x.  = 0 for  all  x 6 X 
(ii)  there  exists  an  x £ X with 

x.  > 0 for  all  j £ N - I. 

J 


Proof ; The  'if"  is  trivially  true  since  (i)  implies  degeneracy. 

Suppose  on  the  contrary  that  (i)  is  not  true,  i.e.,  for  each  j £ N 

there  exists  at  least  one  vector  x^  ^ € X with  xv.^  > 0.  We  can 

J 

choose  a set  of  positive  numbe  -s  X.  such  tr.at  E1?  , X.  =1  and  form 

J j=l  j 

the  linear  combination 


..  I X.x'J> 


J-l 


with  x.  = Z X x^  > 0 
J j=1  J J 
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! 

* 


Therefore  x > 0 and  x £ X .(since  X is  convex),  a contradiction. 

Let  I be  the  maximal  set  satisfying  (i).  I exists  since  it  can  be 

taken  as  the  union  of  all  sets  I satisfying  (i).  If  (ii)  is  not 

satisfied  then  the  set  N - I is  degenerate,  and  by  the  argument 

used  above  there  must  exist  a set  J c N - I such  that  j £ J 

implies  x . - 0 for  all  x£  X,  contrary  to  the  assumption  that  I 
J 

is  maximal.  □ 


( 


( 

< 


The  lemma  shows  that  in  a degenerate  system  there  is  a set 
of  species  which  must  identically  vanish,  while  the  rest  of  the  species 
always  have  a nondegenerate  feasible  composition.  Intuitively,  one 
tends  to  remove  the  vanishing  species  from  the  system  and  consider  a 
canonical  reduced  system  with  species  S^,  i 6 N - I.  The  strong  duality 
theory  of  geometric  programming  justifies  this  reduction  when  F is 
ideal. 

Let  the  reduced  integer  set  be 


r = n - i 


The  problem  formed  by  removing  all  species  which  are  not  in  r is 
called  the  reduced  problem  or  PCPI. 


Problem  rCPI 

Minimize 


i 


Subject  to 


F(x)  = Z x (c  + log  x ) 
j£r  “ “ J 


Z a x . = b . , i £ M 
jer  1J  J 1 

Xj  > o,  j £ r 
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For  convenience  we  shall  use  T as  an  operator,  and  write 
FA  for  the  reduced  matrix, 

I'F  for  the  reduced  objective  function, 
rx  for  the  reduced  composition  vector,  etc. 

Theorem  5.7:  Suppose  problem  CPI  (Section  2.5)  is  feasible  and  degenerate. 

Then  +he  reduced  problem  FCPI  is  canonical.  Furthermore, 

F(CPJ)  * F(FCPI) 

where  F(CPl)  and  F(PCPl)  denote  the  equilibrium  free  energies 
for  the  two  problems,  respectively. 

Proof:  The  boundedness  of  X (Proposition  3.1)  and  Theorems  3 

and  k of  Duff in  et  al.  [25]  establish  the  result.  □ 

■Theorem  5-7  reduces  the  study  of  problem  CPI  to  a study  of 
canonical  programs.  Canonicaiity  can  be  tested  and  the  reduction  per- 
formed using  a linear  programming  method  due  to  Clasen  [13] . His 
technique  finds  a positive  composition  x £ X+  if  one  exists.  It  is 
discussed  further  in  Chapter  7 

A more  difficult  question  is  that  of  degenerate  solutions. 

1 though  an  this  case  one  cannot  perform  an  a priori  reduction,  one 
can  speak  of  a reduced  solution  and  reduced  solution  integer  set  in  a 
way  analogous  to  reduced  systems.  Given  x £ d(f|x),  the  reduced 

solution  integer  se+  is  A(x  ) = l j £ r|x  > 0} . 

J 


lb 


The  A- problem  (A  = A(x  ))  is  the  problem  formed  by  ignoring 

. 

all  species  S.,  j £ N - A(x  ).  We  can  easily  prove 
J 

-ft 

Lemma  $.8:  If  x £ q(f|x)  then 

Ax  £ D(AF|AX) 

where  Ax  , AF,  AX  represent  "reduced"  quantities  and  A = A(x  ). 

* * 

Proof : x £ X implies  Ax  £ AX  since  only  zero  valued 

* 

components  of  x were  deleted.  Suppose  there  exists  a y £ AX 
with 

AF(y)  < AF;Ax*) 

* n 

Define  a vector  y t 1 

yi 

0 

* 

Then  y £ X 

* n * 

F(y  ) = Z y.(c. 

J=i  J J 

= £ * y,(c,  + log  y ) +0 

j€A(x  ) J J J 

= L * y (c  + log  y ) = AF(y) 
j£A( x ) J J J 


j C A(x*) 
j £ N - A(x*) 

+ yt) 

J 
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Similarly 


F(x  ) = AF(Ax  ) 


We  conclude 

ft  ft 

F(y  ) < F(x  ) 

a contradiction,  since  x £ E^f!x)  . □ 

Unfortunately  A is  not  known  in  advance  and  therefore  the 
lemma  cannot  be  applied  to  computations.  Furthermore  A( • ) need  not 
be  unique.  The  preceding  result  is  true  only  for  ideal  F. 

5.4  Uniqueness  of  Solutions 

In  this  section  F is  always  ideal. 

Theorem  5.9:  Let  {F,A,b}  be  canonical,  A £ n > m. 

Let  I =■  {i [there  exists  an  x £ d(f[x)  with  x;.  > 0}  . 

Then  the  solution  to  the  dual  problem  DCP  is  unique  if,  and  only  if, 

rank  (IA)  - m.  Here  IA  = 'A,|j  £ 1}  is  a reduced  matrix. 

J 

Proof ; Let  j I ! denote  the  number  of  elements  in  the  set  I. 

First  we  show  that  if  |l|  = l there  exists  an  x £ d(f|x)  with 
* 

|A(x  );  - l.  The  argument  is  by  taking  convex  combinations,  similar 
to  the  proof  of  Lemma  5-6.  Therefore  we  can  use  A rather  than  I 
and  require  that  there  exist,  an  x £ d(f[x)  such  that  rank(AA)  = m 
where  A = A(x  ) . 
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Theorem  4.3  requires 


* 

u.(t  ) 
J 


* 

that  if  t 
whenever 


solves  DCP  then 

j €A(x*) 


Taking  logarithms  and  writing 


explicitly  we  find 


m 

E 


i=l 


a.  . 
ij 


c 

J 


+ 


log 


j € a(  x 


# 


) 


This  linear  system  in  log  t has  a unique  solution  if,  and  only  if, 

•K* 

the  rank  of  the  coefficient  matrix  is  m.  Note  that  a solution  t 
exists  by  Theorem  4.3  and  that  the  above  linear  system  is  always 
consistent  by  Theorem  4.2.  This  completes  the  proof.  □ 

When  the  conditions  of  the  theorem  are  satisfied,  namely, 
the  solution  to  DCP  is  unique,  the  same  linear  equations  in  log  t 
imply  the  existence  of  a unique  positive  vector  y £ such  that 

A ' log  t - c t log  y . 

Although  y is  strictly  positive,  a nonpositive  solution  x to  CPI 
may  still  exist.  Bigelow  [7]  calls  y tne  vector  of  "virtual  mole 
fractions."  Of  course,  if  |l|  = n,  there  exists  an  x £ d(f|x+), 
i.e.,  x > 0,  and  then  x and  y are  identical. 

As  mentioned  earlier,  the  type  of  uniqueness  implied  by 
Theorem  5-5  does  not  exclude  identical  phases  and  is  therefore 
somewhat  unrealis+ic  for  the  chemist.  As  we  shall  subsequently  show, 
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under  some  regularity  conditions  one  can  obtain  a stronger  uniqueness 
than  was  obtained  in  the  past.  The  regularity  conditions  (on  the  dual 
chemical  problem),  while  mathematically  restrictive,  do  not  pose  much 
difficulty  in  real  systems. 

Consider  the  problem 

Minimize  f(t) 

Subject  to  g^(t)  £ 0,  k = 1,  2,  ...  , K 
where  f(t)  and  g (t)  are  all  differentiable  functions  on  some  set 

A. 

Q c Bn.  This  problem  is  said  to  be  regular  [5]  at  t if 
(i)  gk(t)  < 0,  k = 1,  2,  ...  , K 
(ii)  At  the  point  t,  the  gradient  vectors  Vg^(t)  of  the 
active  constraints  (k  £ K(t)  = (i:g£(t)  = 0))  are 
linearly  independent. 

(iii)  Vf(t)  is  interior  to  the  cone  generated  by  Vgk(t),  k £ K( t ) . 
Conditions  (ii)  and  (iii)  imply  that  all  active  constraints  are 
"strongly  binding,"  that  is,  none  of  them  can  be  relaxed  without 
affecting  the  solution. 

Theorem  5.10;  Let  lF,A,b)  be  canonical  and  F ideal. 

Let  x £ c(F|X)  and  let  A = A(x).  Suppose 
(i)  rank(AA)  = m 

(ii)  The  corresponding  dual  problem  DCP  is  regular  at  an  optimal 
solution  t - t(x).  Then 

x is  the  unique  solution  to  (F,A,b)  . 


Proof:  By  Theorem  5-9>  if  ratik(AA)  = m,  t is  the  unique 


solution  to  DCP.  For  the  dual  chemical  problem,  the  regularity  assump- 
tion means:  T)  > 0 if,  and  only  if,  g (t)  = 1.  Here  tv  , k = 1,2,...,K 
K K K 

are  the  lagrange  multipliers  in  Theorem  h.2  (iii).  From  this  theorem 
we  also  find  that  if  x £ d(F!X) 

xk  = \ = \/&0^  = \ exp[F(x)]  > 0 (1) 

for  all  k £ K(t)  and  x,  = 0 for  k ^ K(t).  Thus,  according  to 

Theorem  $.k,  x.  = 0 for  all  j such  that  k(j)  ^ K(t).  This  is  true 
J 

for  all  x £ D(f|x).  By  Theorem  5 • 5 for  j such  that  k(j)  £ K(t) 
we  have 

x.  = y.  whenever  x,  y £ £)( F | X ) . 

J J 

But  regularity  assures  a unique  vector  T)  of  multipliers.  Thus 
k £ K(t)  are  also  unique  by  (l).  Hence  x = y and  thus,  for  all 

K K 

X,  y s 0(fIx) 

:<j  'VVJ)  = ?A(i)  *yJ  ° 

In  general  we  may  assume  without  loss  of  generality  that 
rank(A)  = m.  We  then  have 

Corollary  5-H:  If  there  exists  an  x £ £)(f[x+)  and  the  dual  problem 

DCP  is  regular  at  a solution  t(x),  then  x is  the  unique  solution  to 

CPI. 
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Corollary  3.12:  If  K = 1 (single  phase)  and  the  system  is  canon- 
ical, it  has  a unique  solution  x = d(f[x+). 

Praof : We  assume  that  rank(A)  = m.  Canonicality  and  bounded- 

ness of  X imply  b / 0 so  that  x £ X implies  x f 0.  In  particular 
x £ D(F'X/  is  nonzero,  so  x > 0 and  by  Theorem  p.L,  x > 0.  The 
result  then  follows  from  Theorem  5-10,  even  without  requiring 
regularity.  □ 


We  have  defined  regularity  in  terns  of  the  dual  chemical  problem. 
It  is  interesting  to  see  what  these  conditions  mean  in  chemical  terms. 
Suppose  that  at  some  x £ £>(F;X)  and  t ^ t(x)  solving  DCP,  the  problem 
is  not  regular.  In  particular,  let  g (t)  = 1 but  q = 0.  Then 

= 0.  Let  e > 0 be  a small  number  and  consider  a change  in  some 
cs,  s £ (i)  such  that 

c = c - e . 
s s 


Therefore  the  constraint  i becomes  binding  and  a solution  of  the 
perturbed  DCP  will  yield  ^(e)  > 0 implying  x (e)  > 0.  Phase  $ 

* a 


\ 

> 


8o 


which  did  not  exist,  appears  after  the  perturbation  for  any  e > 0. 
Hiis  indicates  instability,  or  rather,  discontinuity  in  the  chemical 

system,  reflected  by  discontinuity  in  concentrations. 


As  the  c . 's  are  usually  empirical  numbers  with  finite 
J 

accuracy,  the  probability  of  such  instability  in  practice  is  zero, 
so  chemists  can  be  content  knowing  that  chemical  systems  are  regular 
and  thus  are  likely  to  have  a unique  equilibrium  composition. 


5.5  Extensions  to  Nonideal  Systems 

Geometric  programming  duality  serves  well  in  analyzing  the 
ideal  case.  Unfortunately,  no  parallel  theory  was  developed  for  the 
more  general  nonideal  case.  Moreover,  there  is  no  guarantee  that  a 
meaningful  duality  even  exists.  In  this  section  we  shall  assume  only 
homogeneity  of  degree  one  and  convexity,  as  postulated  in  Chapter  3. 

We  shall  examine  some  conditions  leading  to  dual  problems,  particularly 

»»  II  _ - 

pure  duals . 

A major  complication  in  subsequent  analysis  is  that  F(x) 
may  be  undefined  and  surely  is  not  differentiable  when  x ^ 0.  We 
shall  bypass  these  obstacles  by  sacrificing  some  generality  and  dealing 
only  with  positive  compositions.  Our  basic  problem  is  CPN  (Section  2.5). 


Theorem  5.15:  Let  F be  differentiable  on  1^,  convex  and  homo- 

geneous of  degree  one.  Consider  problem  CPN.  Ihen  there  exists  an 

x Z £(f|X+)  if,  and  only  if,  there  exists  a vector  z £ tT  su  h 
* 

that  z solves  the  problem « 
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Maximize  z'b 

•* 

Subject  to  z'A  <(i'  (x 
Furthermore, 


* . * * * 
z b = F(x  ) and  A z = p(x  ) 


Here 

F(x)  = |r’(x)-x 

n(x)  is  the  vector  of  chemical  potentials. 


Proof;  By  direct  application  of  Kuhn-Tucker  conditions  [36] 
to  the  convex,  lineariy  constrained  problem  CPN,  the  necessary  and 
sufficient  conditions  based  on  the  Lagrangean 

L(x,z)  = F(x)  - z'(Ax  - b) 

are 

(i ) Ax  = b 

(ii)  V^L^x  >z  ) > 0,  i.e.,  p(x  ) - z A > 0 
. . * * , *,  * * 

(iii)  z Ax  = p(x  j*x  = F(x  ) (complementary  slackness). 


From  (i),  (ii) 

* * * * * * 
p(x  )-x  = F(x  ) > z Ax  = z b 

From  (i),  (iii)  and  the  condition  x > 0 

* * * . *.  , * . * 
z Ax  = z b - F(,x  ),  A z = p(x  ) 


Tne  rest  is  clear.  □ 
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The  dual  problem  implied  by  the  theorem  is  not  "pure"  for  it 
involves  primal  variables  x.  Eisenberg  [28]  studied  the  existence 
of  pure  dual  problems  to  homogeneous  programs"  (■where  F(x)  is 

homogeneous  of  degree  one),  but  his  assumptions  are  not  valid  in  our 
case . 


Let 

F(x)  = p'(x)’x  and  suppose 

(1) 

“j(x)  1 

= Hj(Xj),  j £ (k),  k = 1,  2,  ...  , K 

(2) 

-j<J) 

is  well  defined  for  x.  > 0,  j £ N 

t} 

(3) 

t;V 

is  strictly  monotone  in  * and  p'^.y)  exists  for 

0 J 

all  y 

in  the  row  space  of  A. 

Then  the  constraints  of  Theorem  5.13,  namely, 

A'z  < u(x*) 

can  be  inverted 

u.^A  !z)  < x* 

J J ~ j 

Summation  over  j £ (k)  yields 

Sk(z)  = Z »x"1(A  ’z)  < 1 
J£(k)  J 

Therefore,  if  ,(l),  (2)  and  (3)  hold,  a pure  dual  resembling  a geometri 
program  can  be  formed.  The  hardest  condition  to  accept  is  probably 
condition  (3).  Condition  (l),  although  quite  restrictive,  is  never- 
theless common  in  ideal  or  nearly  ideal  systems.  It,  is  natural  to 
expect  that  p(x)  = p(x)  since  p is  homogeneous  of  degree  zero. 


In  practical  nonideal  systems,  p(x)  is  based  in  many  cases 
on  empirical  relations  (see  Section  8.7  for  example)  which  are  valid 
in  some  range  of  x but  need  not  preserve  either  homogeneity  or 
convexity.  The  best  one  can  hope  to  do  in  the  general  situation  is 
to  be  able  to  handle  such  problems  computationally. 

We  conclude  our  short  treatment  of  the  general  case  with 

Proposition  5-14 : Let.  F(x)  = n(x)  •x  be  differentiable  on  E^, 

convex  and  homogeneous  of  degree  one.  Consider  problem  CPN.  Then 

x £ fi(F | X+) 


if  and  only  if 


e'n(x)  = o 


for  every  reaction  vector  8. 

Proof:  By  Theorem  5- 13  there  is  a z such  that 

z 'A  - p(x) 


Consequently 


z 'A 9 = p(x) "9 


but  by  definition  of  reaction  vectors  At?  = 0.  Thus 

0 = p(x)-<?  □ 
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The  relation  n(x)-S  = 0 is  the  analog  to  the  logarithm  of 
the  mass  action  laws  (Section  3-2).  Restoring  the  usual  form,  we 
can  write  the  generalized  mass  action  laws  as 


In  the  ideal  case: 


n 6 u (x) 
n e J 3 = 1 

j=l 


li.^x)  = c,  + log  x . We  have 
J J J 


n 

n 

J=i 


exp(-c ' •&)  - K(e) 


as  was  shown  in  Chapter  3- 


CHAPTER  6 


APPLICATIONS  OF  CHEMICAL  DUALITY 

Several  applications  of  the  duality  presented  in  Chapter  4 will 
be  demonstrated  in  this  chapter.  A major  application,  computational 
algorithms  based  on  duality,  will  be  discussed  in  the  next  chapter. 

The  chemical  system  is  {B.S,#},  characterized  mathematically  by 
(F,A,b)  as  described  in  Chapter  2.  The  applications  presented  are 

a.  Testing  for  equilibrium,  given  a feasible  composition. 

b.  Verification  of  the  model- -upper  bounds  on  concentrations  of 
trace  components. 

c.  Sensitivity  ana lysi s- -changes  of  the  solution  with  variation 
in  the  free  energy  parameters. 

d.  Goaling  techniques-- solution  of  problems  with  side  conditions 
("goals")  on  the  equilibrium  composition. 


6.1  Testing  and  Characterization  of  Equilibrium 

Chemical  duality  can  be  applied  in  n simple  and  straight- 
forward way  to  test  whether  a given  composition -vector  x £ X is 
an  equilibrium  solution  for  an  ideal  system.  Eie  results  of  the 
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previous  Chapter,  especially  Theorem  5*4,  indicate  that  if  x £ d(f|x) 
(the  "equilibrium  set")  then  x^  = 0 implies  = 0,  so  we  assume 
that  this  condition  is  satisfied.  For  an  ideal  system,  geometric  pro- 
gramming duality  Theorem  4.2  and  the  form  of  the  dual  chemical  problem 
DCP  indicate  that  x £ d(f|x)  if  and  only  if  there  exists  z £ ^ 


such  that 


a,  .z.  = c . + log  x.A  for  j £ N such  that 
ij  1 J j 


Vj) 


> 0 (l) 


This  equation  can  be  easily  obtained  by  applying  the  Kuhn-Tucker  Con- 
ditions (see  [7])-  Thus,  an  optimality  test  will  consists  of  solving 
the  linear  system  (l)  for  z.  In  any  nontrivial  case,  the  number  of 
equations  willexceed  the  number  of  variables,  so  that  this  is  a test 
for  the  consistency  of  system  (l).  A measure  of  how  far  x is  from 
the  optimum  can  be  obtained  by  finding  the  "least  squares"  solution 
to  the  system.  Let 

A = A(x)  = (j  | x . > 0) 

J 

AA  = (A  . | j £ A) 

J 

Ax  = | j £ A) 

Ac  = (c j | j £ A) 

Then  the  least  squares  solution  z is  given  by 

z = IaA(AA)']  "''•AA  •{log(Ax)  + Ac] 

whe~°  AA  is  assumed  to  be  of  full  rank  m.  If  x is  an  equilibrium 


solution  then  the  error  e is 


e = (AA)'z  - log ( Ax)  - Ac  = 0 
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Otherwise,  e is  an  appropriate  measure  of  deviation  from  equilibrium. 
2 

Note  that  e is  actually  a measure  of  the  deviation  from  a solution 
which  satisfies  the  mass  action  laws. 


b .2  Verification  of  the  Model 

This  application  was  first  noted  by  White  et  al.  [5^].  We 
expand  their  result  by  adding  bounds  on  concentrations  of  added 
species.  A model  (E,S,G>}  is  based  on  the  hypothesis  that  S is 
indeed  the  right  set  of  species  and  that  all  other  possible  species 
may  appear  only  in  negligible  amounts.  The  purpose  here  is  to  test 
this  hypothesis. 

It  is  assumed  that  an  equilibrium  composition  x for  the 
model  is  known,  together  with  the  dual  chemical  variables  z computed 
by  (l.l).  We  wish  to  test  whether  the  exclusion  of  some  species 
^ S is  justified.  We  assume  that.  k(q)  £ $ and  "=*  0. 

Let  A = la.  ) be  the  formula -vector  of  S . By  (l.l)  we  must  have 

q iq  q 

at  equilibrium 


m 

z 


i=l 


a . z . 

iq  l 


c 

1 


(1) 


Thus,  the  approximate  concentration  of  S , had  it  been  included  in 

Q. 

the  model,  would  be  x' 

q 


xJ  = exp[  Z a.  z. 
q i.,1  iq  1 


Cq] 


(?) 
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If  x®  is  small  (say  x^  < 10  ^),  one  can  usually  justify  the  exclusion 

q Q. 

of  S • The  true  value  of  x , had  S been  included,  is  of  course 

q q q 

different  since  with  S in  the  model,  z is  no  longer  optimal.  In 

q 

many  situations,  one  can  obtain  an  upper  bound  on  any  x^  (whether 

S € S or  not)  without  ever  solving  the  problem,  by  noting  that  since 
J 

x,  < 1,  we  always  have 
J 


,2,  auzi  £ V 


j * 1,  2, 


when  z solves  DCP.  To  find  a bound  on  x^  we  solve  the  linear  pro- 


Maximize  y = E a.  z. 

q i=1  iq  i 


Subject  to  A'z  < c 


Letting  y be  the  maximum,  we  obtain  by  (2) 


S exP[yq  - cql 


The  right  hand  side  may  be  greater  than  unity,  in  which  case  the  bound 

* 

is  useless.  This  happens  when  y > c . Solving  a linear  program  to 

q q 

obtain  a single  bound  is  of  course  highly  inefficient.  In  many 
practical  situations,  upper  bounds  can  be  found  by  inspection.  Recall 
that  there  is  no  loss  of  generality  in  assuming  that  A > 0 for  a 
real  system  (ignoring  electroneutrality),  and  that  each  subspecies 
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is  also  included  as  a species  where  j(i)  is  the  index  of  the 

"species"  • Using  this  idea  in  (3)  we  have 


Zi1Cj(D 


(5) 


because 


1 I = i 


0 otherwise 


If  A > 0 then  for  any  J,  and  any  z solving  DCP 


m m 

E aijZi  = A aiJCj(i) 

i=l 


(6) 


Similar  to  equation  (4),  we  arrive  at 


xq  < exp{yq  - cq] 


<7) 


Again,  if  y^  > c^,  the  bound  is  useless  since  the  right  hand  side  is 
greater  than  one.  However,  we  found  this  technique  very  useful  in 
obtaining  bounds  and  initial  approximations  with  the  dual  algorithm  of 
the  next  chapter. 


Example  6.1.  Suppose  one  wishes  to  test  whether  ozone  0^  is  justly 
excluded  from  the  hydrazine  model  (Appendix  A. 3).  It  is  assumed  that 
this  could  be  species  S^  with  ~ -14  (approximated  for  3500° 
from  data  in  [38]).  The  formula  vector  for  0^  is  (0,3>0),  corre- 
sponding to  subspecies  (H,0,N).  Now  observe  that  Sq  = oxygen  (0) 
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and  cq  = -10.708.  Thus,  zg  < -10.708  from  equation  (5).  Then, 
by  (6) 


yll  = 3'c8  = -32-124 


Equation  (7)  leads  to 


xn  < exp[-32.124  - (-14)] 

^ . -18.124  , -7 

$11  < e <10 


The  small  concentration  indicates  that  0^  is  only  a trace  element 
and  can  be  excluded  without  affecting  the  model. 

Up  until  now  it  was  assumed  that  the  phase  of  the 

excluded  species  S^,  exists  at  equilibrium.  Difficulties  arise  when 


the  phase  does  not  exist,  i.e., 


*k(q) 


= 0 at  equilibrium.  We  wish 


to  test  whether  inclusion  of  S in  the  model  would  have  caused  the 

q 

phase  to  appear. 

Assume  that  the  dual  chemical  problem  DCP  is  regular  at  a 
dual  solution  point  z (Section  ^.b).  For  the  dual  chemical  problem 
DCP  (without  S ) we  have 


0 

Sw (a)  = £ exPf  £ a. ,z  - c ] < 1 . 

J€(k(q))  i=l  iJ  1 j 


After  inclusion  of  S we  obtain 

q 
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(8) 


B*U>  ■ 4(1)  + Vi  ’ V 

= 0 and  the  phase  still  vanishes. 

Otherwise,  the  phase  appears,  in  which  case  the  concentrations  of  species 
3 j > Je*k(q)»  can  be  approximated  by  Uj(z)/gj^qj,  i.e., 

m 

= — - — • exp[  Z a z - cl,  j € (k(q)> 

J «k(q)  1=1  J J 

Ihe  vector  z is  no  longer  the  exact  optimal  solution  to  the  dual 
DCP  of  the  expanded  problem  with  If  is  no't  small,  the 

mass  balance  equations  are  violated  and  one  has  to  resolve  the 
problem  with  S^  included  in  the  model. 


if 


e*(q) 


<.  x tnen 


^(q) 


6.J  Sensitivity  Analysis 

When  the  parameters  of  a given  system  are  changed,  its 

equilibrium  composition  may  also  change.  Bigelow  and  Shapiro  [9] 

analyzed  the  variation  of  the  equilibrium  composition  x with  small 

changes  of  A,  b,  and  c by  using  geometric  programming  arguments. 

A limited  earlier  version  appeared  in  Duffin  et  al.  [25] . Later 

analysis  by  Theil  [49]  was  concerned  with  variations  in  primal 

geometric  variables,  the  variables  of  the  dual  chemical  problem. 

For  most  practical  purposes,  the  changes  of  interest  occur 

primarily  in  the  free  energy  coefficients  c.  when  the  system  is 

J 

ideal.  A special  interesting  case  occurs  when  the  vector  c changes 
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with  a single  parameter,  the  temperature  T for  example.  In  bio- 
chemical systems,  where  changes  are  usually  small  due  to  internal 
stabilizing  mechanisms,  sensitivity  analysis  can  analyze  corresponding 
changes  in  composition  without  solving  the  problem  afresh. 


Computing  (l),  we  arrive  at 


K 

-b.  + E T).[  E a u (z)]  =0,  i = 1,  2,  ...  , m (5) 

1 k=l  k J€(k)  J 


Let  c + ac  be  the  new  coefficient  vector  with  c being  a direction 

vector  and  a a parameter.  The  corresponding  changes  in  tj,  u(z), 
»*/>*!»* 

and  z are  ocrj,  Ou,  az,  respectively,  a can  be  ignored  since  our 
interest  is  only  in  directions  of  change.  We  assume  that  these 
directional  derivatives  exist,  at  the  solution  point  z.  (Existence 
requires  that  the  conditions  of  the  implicit  function  theorem  be 
satisfied  for  problem  DCP.  Biis  topic  is  thoroughly  studied  in  [9] . ) 
Conditions  (l)-(4)  must  hold  also  for  c + ac,  tj  + ar),  u + au, 
z + az,  so  that  (5)  leads  to 


E T]-[  E a u,(z)]  + E T).  •[  E a u (z)]  = 0 


jeu> 


iJJ 


je(k) 


ij  T 


(6) 


From  the  definition  of  u (which  is  considered  a constraint) 


“j(z)  s tfaiA1'uj<2) ' Yj(z) 


(7) 


If  g^(z)  = 1,  we  can  write  by  (2) 


From  (3) 


(9) 


\*°  if  \ = 0 

and 

\’(i  - g^U))  - \‘ek(z)  = 0 (10) 

Assume  that  the  problem  DCP  is  regular  at  z.  Ihen 


= 0 implies  g^z)  < 1 

In  this  case,  if  g (z)  = 1 then  (10)  implies  gfc(z)  = 0 so  that 
(8)  is  redundant.  Similarly,  if  tj  (z)  = 0 then  (10)  implies 
t^(z)  = 0 so  that  (9)  is  redundant. 

For  the  regular  case,  therefore,  (6),  (7),  and  (10)  cover 
all  the  conditions. 

Let 


and  let 


hik" 


z 

j£(k) 


auuj(z) 


Combining  (6)  and  (7)>  changing  order  of  summation,  and  inserting 
the  definitions  of  h and  d leads,  after  some  algebraic  manipu- 
lations, to 


+ Z a 

j 


ijdj 


(Z 


ij 


a d c 

ij  J j 


(11) 


for  i = 1,  2,  ...  , m 
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Equation  (10)  leads  to 


J€(k) 


•k  "J 


u,(z) 


\ ? hik'zi  + 


je<k) 


Vj  - 0 


for  k = 1,  2, 


, K 


Let 

A = (aij},  ASiT1 
D = diag(dL,  dg,  ...  ,dn>- 
H = (hik),  H £ 

»V  )V  K 

n = (nkL  r\e  tf 

A = diag(Tj1,  • • • ) T)jj) 

- (1,  1,  , 1)  € / 


G(z)  = diagonal  matrix  in 


■P'S  0^2)  = 1 - g^U) 


E = (ekJ},  E € IKXn 

and 


j € (k) 
otherwise. 


With  these  definitions  and  some  additional  algebra  we  arrive  at 


-AH  i G(z) 


This  is  a linear  system  of  m+K  equations  in  m+K  unknowns— z and 
i).  If  this  system  is  nonsingular  (which  we  assume,  of  course),  z 
and  T)  can  be  computed  explicitly  for  any  giver.  (provided  that 
the  solution  at  c,  namely  z and  t),  is  known).  To  recover  the 
composition -vector  x,  our  main  interest,  recall  that  Uj(z)  = x^ 
whenever  > Applying  equation  ( 5 ) we  obtain  the  mass 

balance  equations 


z 

j 


i = 1,  2, 


m 


whence 

Vnk(j)  = Xj  = Xj‘*k(j) 


Therefore 


uj’\(j)  + Wj) 


(1*0 


(15) 


u,  can  be  computed  from  z via.\7). 
J 

(15). 


T)  is  found  directly  from 


In  conclusion,  we  have  found  the  directional  change  x due 
to  a change  c in  the  coefficients. 

We  shall  now  examine  a special  case  where  c changes  with  the 
temperature,  i.e.,  c = c(T)  and 


~ _ dc(T) 

dT 


Recall  that  by  equations  (2.4.5),  (2.4.6)  in  the  ideal  case 
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(16) 


ti?(P,T) 


(P,T)  = ^ RT  + log[r,(P,T)] 


We  assume  that  r.(P,l)  is  independent  of  T and  that  P is  fixed 
J 

so  that  r (P,T)  - r (constant).  From  elementary  thermodynamics 
J J 

[21,  p.  lkj>] 


d[p1(x,P,,T')/RT] 


dT 


AE.(T) 

RT‘2 


(17) 


Here  AH,(T)  is  the  enthalpy  of  formation  of  S at  T.  R is 
J J 

the  gas  constant.  Furthermore 


d[£H  ^(T)} 
dT 


* cpj(t)  : qj  + 6T  + rf 


(18) 


where  C ^ is  the  heat  capacity,  a,  P,  and  r are  empirical  coefficients. 

We  can  integrate  the  last  relation,  adding  an  integration  constant  H . 

J 

Substituting  in  (17),  integrating  again,  and  identifying  ^ with 

J 

c.  (for  the  pure  species),  we  obtain 
J 


1 0 

c , - c . 
j J 


tr)  - a -iog(-i) 

io  J o 


_ _i  • ' ^ - T2)-  — 1 ‘(T'>  - T^)l 
4 x 1 01  9 U1  CrJ 


where  the  integration  is  from  initial  aboslute  temperature  to 

final  temperature  T.  and  c^?  - c (T  ).  For  small  changes  of  T, 

r J J ^ 
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like  those  expected  in  biological  systems,  the  first  term  is  usually 
dominant.  Let  H = (H-^,  Hg,  ...  > Hn)  and  taking  the  < finite) 
difference  c^  - c^  as  direction  we  obtain  the  well  known  formula 


~ H 


_ ,_1_  ±s 

•»  ' m “ T 


R T, 


Values  of  H are  typically  tabulated  in  handbooks  of  chemistry  for 
the  temperature  T = 298.l6°K.  Changes  in  temperature  can  thus  be 
directly  translated  to  changes  in  composition. 


6.4  Goaling  Problems 

Goaling  problems  are  variations  of  the  classical  chemical 
equilibrium  problem.  They  differ  from  the  standard  problem  in  that 
some  predetermined  conditions  ("goals")  are  imposed  on  the  equilibrium 
composition.  An  attempt  to  solve  the  problem  by  simply  adding  the 
goals  to  the  mass  balance  constraints  will  violate  the  principle  of 
internally  unconstrained  equilibrium,  discussed  in  Section  3.5.  To 
compensate  for  the  loss  of  internal  "degrees  of  freedom”  due  to  the 
goals,  one  normally  relaxes  some  of  the  "external"  mass  balance  con- 
straints. This  is  done  by  allowing  an  open  system,  in  which  some 
species  can  be  freely  added  externally  to  achieve  an  equilibrium 
where  the  goals  are  satisfied.  Although  goaling  problems  do  not 
resemble  our  familiar  chemical  equilibrium  problem,  duality  can  i 1 
some  instances  modify  the  problem  to  put  it  in  a standard  chemical 

99 


m 


ijimui,)iiiw«i^ 1 WWW  I u.itnjy. ' jy « 1'  W.WIV  .»Trri  '*'*  •"•*• 


equilibrium  format.  Our  main  result  here  is  a generalization  of  a 
procedure  for  concentration-goaling  by  DeHaven  [19]  . Ihe  problem 
is: 

Problem  CPG 

Given  an  ideal  rystem  (F,A,b),  find  the  equilibrium  composition 

and  the  amount  of  species  to  be  added  (or  removed)  externally, 

such  that  the  equilibrium  concentration  of  Sr  is  some  given  value 

DeHaven  [19]  and  the  RAND  chemical  composition  code  [48]  deal 

only  with  the  case  where  the  species  appears  also  as  a subspecies. 

We  shall  relax  this  assumption.  Furthermore,  we  show  how  to  extend 

the  procedure  to  several  simultaneously  goaled  species.  It  iB 

assumed  that  the  problem  is  canonical  and  has  a solution.  Hie  formula- 

vector  of  S is  A . 

r r 


Theorem  6,2;  Problem  CPG  is  equivalent  to  the  chemical  equilibrium 
problem  CP#  defined  by  the  triplet  {F#,  A#,  b#}  where 


F#  - F + [-(cr  + log  xr)  + log  xn+1]xn+1 


A#=(a#}  = 


b#  = (\)f)  = b - 0 

1 i 


a,  , 

1 < 

u 

”air 

J = 

- 9 • a 

ir 

Kr 

> 0) 
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Proof : First.,  note  that  CP"  was  constructed  from  CPG  by 

adding  a species  Sn+1  which  is  in  a separate  phase  ®K+1“  CP^ 
is  canonical  since,  by  hypothesis,  the  underlying  ungoaled  program 
CP  is.  To  see  this  let  x £ X+(A,b).  For  any  e > 0,  let 


/ > 0 


xr  + e j - r 


* < s ♦ 9 j = n+1 


otherwise 


A#.x#  = Ax  + A e - A *(e  + 6) 

■yi  k ' 


r r 


= Ax  - A fi  = b - A 6 = b# 
r r 


/ € X#V,b#; 


Next,  observe  that  if  i - arg  min{'bi/air Iair  > C],  then  x^/a^ 

equals  the  total  number  of  moles  of  subpsecies  l,  which  is  unknown 

due  to  the  open  ended  nature  of  the  problem. 

* * 

If  problem  CPG  has  a solution  x by  adding  xn+1  moles  of 

* * 

Sr,  then  there  is  a solution  to  problem  CP  defined  by  (F,A,b  } 

* -fc- 

where  b = b + A *x  . The  solution  x satisfies  the  goal: 
r n+l  ° 

.V*  A 
X = X . 

r r 


10.1 


Suppose  x were  known,  then  by  duality,  there  exists  a 
n+i 

vector  z such  that 


fAj  = Gj  + log(xy  fcr  a11  J such  that  > 


In  particular 


z 'A  - c + log(  x 't 
r r > r 


* * 

From  the  solution  x to  CP  we  cons+ 


truct  a solution  x $ to  CP^ 


as  follows 


4 - < 


1 < j < n 


0 + x j = n+1 

n+1 


Indeed 


J.v# 


AF  ■ x77"  = Ax  - A ' ( x id) 
r n+i 


-*  M. 

= b + A x , • A 0 - A x br 

r n+^  r r n+x 


# * 

Clearly,  if  x solves  CP  then  z above  will  satisfy 

-z'A  = -(c  + log  x ) -r  log  1.  Therefore  $ is  a solution  to 
r r r 

CF^.  Conversely,  any  solution  to  CP^  must  have  a dual  chemical 

solution  z such  that 


•z'A  = -( c +•  log(x/ ) 
r r r 


Hence 
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log  = -c  + z'A  = log  X 
r r r r 

It  is  clear  from  the  construction  that  x^  will  satisfy  the  mass 
balance  constraints  and  x^+^  ^ 0.  If  x^f+^  = 0 all  species  with 
a.  . > 0 must  vanish,  but  we  can  assume  without  loss  of  generality 

1.  Jj 

that  at  least  one  such  species  exists  in  each  phase.  Theorem  5.4 
then  leads  to  x.  = 0 for  all  j,  which  is  contrary  to  the  hypothesis 

J 

that  a solution  exists.  This  completes  the  proof.  □ 

Corollary  6.3:  The  procedure  indicated  by  the  theorem  can  be  extended 

to  simultaneous  goaling  of  several  species. 

Proof:  One  can  proceed  inductively,  goaling  on  species  S 

1 r 

generates  problem  CP^.  Starting  with  CP^  and  goaling  on  S 

s 

generates  (CP#)#  and  so  on.  □ 

At  most  m independent  species  can  be  goaled  this  way,  but  even  a 

smaller  number  may  lead  to  inconsistencies. 

We  are  unable  to  find  a similar  method  for  the  more  general 

problem  of  goaling  S by  adding  or  removing  a different  species  S . 

I*  s 

It  is  possible,  to  cast  this  problem  in  the  form  of  a parametric 
chemical  equilibrium  problem  with  a single  parameter  a. 

Let 

£ 5 arg  min(b./a  |a  >0} 

1 IS  IS 

and  let 

8 S Va.ls 

10} 


r— v., 


»-..v 


p 


s 

¥ 
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j 

i 


I 

■( 

f 


& is  interpreted  as  the  total  amount  of  Sg  currently  in  the  input. 

* 

Set  b = b - 9A  where  A is  the  formula  vector  of  S . Let 
s s s 

a > 0 be  a parameter  and  set 

b^=  b^(ct)  = b*  + aAg 

The  problem  is  now  characterized  by  (F,A,b$(aO) . To  solve  it, 

one  has  to  find  a,  and  the  equilibrium  solution  x(a),  such  that 

x (a)  will  have  :he  foaled  value.  The  amount  of  S to  be  added 
r s 

(or  removed)  is  ol-Q , assuming  that  9 is  the  original  input  of  S . 

s 

The  class  of  goaling  problems  can  be  generalized  further  by 
requiring  the  equilibrium  composition  x to  minimize  a given  function 
cp(x).  Again,  this  goal  is  achieved  by  allowing  an  open  system  with 
varying  input  The  problem  can  thus  be  formulated  as  follows 

Problem  CP(fl) 

Let  fi  be  a subset  of  if1. 

Find  y £ Q and  a composition-vector  x = x (y) 

x*  € X(A,y) 

such  that 

x*  £ D(F|x(A,y)) 

and  cp(x  ) < cp(x)  for  all  and  x,  such  that 

x € d(F|X(A,  w)). 

As  a simple  example,  consider  the  case  where 
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fi  = (y  e #|y  > 0,  Hyll  = 1) 

and 

cp(x)  = x"1  . 

In  chemical  terms,  we  wish  to  find  the  input  vector  y (replacing 

the  usual  fixed  b)  such  that  x , the  concentration  of  S , is 

r r 

maximized.  The  solution  is  trivial--an  input  of  pure  Sr»  To  see 
this,  notice  that  in  the  dual  chemical  problem  DCP(fl)  associated 
with  CP(fl),  nothing  changes  except  the  objective  function  which 
now  reads 

minimize  -y'z  . 

x is  identified  with  the  term  u (z) 
r r' 

u (z)  = exp[Z  a z - cl 
r ^ ir  i r 

To  maximize  xf  we  can  maximize  u^(z)  in  the  dual  chemical  problem. 
The  problem  is  then:  find  y 6 ft  such  that  the  solution  z(y) 

to  DCP(n)  will  also  be  a solution  to  max^Eu^Ez))  where  z = z(y). 

The  goaling  condition  (in  the  dual)  can  be  written 

minimize  t(z) 

* 

By  choosing  y^^ 


= - Z 

i 


a z 
ir  i 


ir 
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one  can  easily  verify  that  y £ fl  and  that  any  solution  z to  the 

* * 
dual  problem  of  CP(y  ) will  also  minimize  \jr(z).  In  this  case  y 

was  found  even  without  solving  CP(fl).  The  equilibrium  composition 

can  now  be  found  by  solving  the  problem  CF(y  ),  characterized  by 

(F,A,y*} • 

A useful  example  of  this  type  is  when 

cp(x)  = X /x 
r s 

Equivalently,  we  wish  to  maximize  x /x  . This  particular  application 

is  of  interest  when  the  equilibrium  products  have  to  be  separated, 

for  instance,  when  the  solution  yields  some  undesired  byproduct. 

The  cost  of  separation  depends  significantly  on  the  concentration 

ratio  of  the  species  involved.  Naturally,  one  wishes  to  maximize 

the  ratio  of  desired  product  to  undesired  product. 

In  an  analogous  way  to  the  method  described  before,  the 

problem  is:  minimize  y(z)  = 2. (a.  - aJ  )z.  over  all  y £ ft  such 

i l r l s l 

that  z solves  the  dual  problem  DCP(y),  i e.  , over  all  y £ fj  such 
that  z solves 

Minimize  -z'y 

Subject  to  g^(z)  < k » 1,  2,  ...  , K. 

Unfortunately,  here  The  solution  is  not  immediate,  thus  some 
iterative  procedure  is  needed.  Goaling  problems  are  closely  related 
to  a problem  treated  by  Dantzig  et  al.[l6]  and  to  the  so-called 
"pooling  problem." 
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In  general  one  can  apply  the  Kuhn-Tucker  conditions  to  create 
a mixed  primal-dual  problem, which  for  problem  CP(fi ) with  goal  function 
cp(x)  is 

Minimize  cp(x) 

Subject  to  kx  - y 
x > 0 
y £ a 

A 1 z = c + log  x 

This  form  is  not  very  attractive,  as  it  resembles  neither  the  primal 
nor  the  dual  problems  Most  of  the  useful  structure  of  either  primal 
or  dual  problems  is  lost.  Finally  we  note  that  with  goaling,  the 
assumption  of  an  open  system  implies  that  X is  no  longer  bounded. 
Most  of  our  previous  results,  especially  those  which  assume  existence 
of  a finite  x ££)(f!x)  are  still  in  force  for  the  unbounded  case. 
Unbounded  problems  were  treated  by  Bigelow  et  al.  [8],  and  by 
Kortanek  et  al . [ 3 5 J - 
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CHAPTER  7 


A DUAL  CHEMICAL  ALGORITHM 


7.1  Introduction  and  Review  of  Existing  Methods 

Most  work  on  tne  chemical  equilibrium  problem  has  been  directed 
toward  development  of  computational  algorithms.  With  increasing 
dimensions  and  complexity  of  the  problems  treated,  neither  ad  hoc 
techniques  for  specific  problems,  nor  the  time-honored  approach  of 
solving  a linear-nonlinear  system  of  mass  balance  and  mass  action 
equations  are  adequate  any  more.  Van  Zeggeren  and  Storey  [53]  present 
a comprehensive  review  of  techniques,  divided  into  two  broad  classes, 
one  using  the  classical  approach  of  solving  the  mass  action  equations, 
the  other  using  optimization  techniques.  Notable  in  the  first  class 
is  the  pioneering  work  of  Brinkley  [10]. 

Among  the  optimiza+ion  techniques  we  mention  Newton-type  methods 
based  on  quadratic  approximations  by  White  et  al.  [5^]  and  Dluzniewski 
et  al.  [22],  the  separable  programming  approach  of  Dantzig  et  al.  [l8], 
Dantzig's  generalized  programming  method  [15],  Clasen's  linear- 
logarithmic  algorithms  [12-lL],  and  Bigelow's  linear,  quadratic, 
and  dual  algorithms  [7] • All  these  methods  can  be  classified  as 
primal  chemical  algorithms  since  they  all  work  directly  with  the 


108 


chemical  problem,  improving  the  composition  vector  at  each  i ration. 
Some  of  the  primal  chemical  algorithms, especially  Dantzig's  [15], 
Clasen's  [12-lb]  and  Bigelow's  methods  [7]  use  dual  multipliers  within 
each  iteration,  but  are  nevertheless  primal  oriented.  Passy  and 
Wilde  [bl]  were  the  first  to  implement  a pure  dual  chemical  method, 
but  their  approach  is  valid  only  for  single-phase  problems 
Bigelow  presented  two  "dual"  metnods  and  suggested  some  of  their 
advantages,  without  reporting  any  computational  experience. 

After  the  relation  of  chemical  equilibrium  to  geometric  pro- 
gramming was  realized,  geometric  programming  problems  were  solved 
using  chemical  equilibrium  codes,  especially  the  highly  efficient 
program  developed  at  RAND  [b8]  In  a sense,  the  algorithm  presented 
here  does  just  the  opposite--it  uses  a geometric  programming  method 
to  solve  chemical  problems . 

Two  points  motivated  the  study  of  this  dual  chemical  algorithm: 
(i)  The  dimensions  of  the  problem  are  significantly  reduced,  from 

n variables  and  m constraints  in  the  primal  chemical  problem, 
to  m variables  and  K (nonlinear)  constraints  in  the  dual. 
This  reduction  can  be  advantageous  both  in  computing  time, 
and  certainly  in  the  storage  requirements . 

(ii)  The  inability  of  the  RAND  code,  and  of  most  other  primal 
chemical  methods  to  handle  degeneracies  in  the  sense  of 
Chapter  5.  Problems  with  vanishing  phases  and  linear  de- 
pendencies due  to  ill-formulated  problems  are  the  two  main 
causes  of  such  degeneracies. 
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As  our  basic  method,  we  have  chosen  an  algorithm  suggested 
by  R.  Dembo  [20] , originally  developed  to  solve  geometric  programming 
problems  with  both  positive  and  negative  coefficients.  Since  the 
dual  chemical  problem  involves  only  posynomials,  we  shall  be  inteiested 
only  in  the  standard  geometric  programming  part  of  his  technique.  In 
its  original  form  the  algorithm  was  found  unsuitable  for  chemical 
problems  due  to  the  large  coefficient  values  in  the  dual.  After  pre- 
senting the  basic  features  of  the  algorithm  we  shall  present  a modified 
version  of  it,  shown  to  be  equivalent  to  Zangwill's  concave  cutting 
plane  method  [55]-  The  theory  is  followed  by  some  computational 
results  and  a discussion  of  computational  aspects. 


7.2  The  Geometric  Inequality  and  Condensation 

A general  form  of  the  well  known  inequality  between  the 
arithmetic  and  geometric  weighted  means  of  positive  numbers  states 
that 

Theorem  7-1:  Let  u.. , u_,  ...  , u be  positive  numbers  and  let 

1 d n 

w, , w^,  ...  , w be  nonnegative  numbers  such  that  Z?  , w,  = 1. 

1 27  n i=l  i 

Then 

n n w 

£ w u.  > II  u. 1 (l) 

i=l  1 " i=l  1 

Furthermore,  equality  holds  if,  and  only  if,  there  is  a constant  K 
such  that  u^  = K for  all  i.  A proof  can  be  found  in  Hardy,  Little- 
wood,  and  Polya  [ • 
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This  inequality  called  henceforth  the  "geometric  inequality" 
plays  a key  role  in  geometric  programming  It  will  now  be  put  in 


a more  general  form,  where  the  u^  s and  w_^'s  are  positive  and 
nonnegative  functions  respectively. 


Theorem 


7.2.  Let  B * (x  £ Bjx  > 0},  R lx  £ l| x > 0} . Let 


JL 

u^n^  c JT  — > 1+,  w^  ‘fig  c l1  — > fi  , 1 = 1,2,  . . , , n,  be  functions, 

such  that  Z1}  , w.(x)  = 1 for  some  x £ f?o. 

1=1  i c / \ 

w (x) 

When  w (x)  = C we  define  [w  (x)'J  1 = 1.  Then 


n n pu.(t)-jWivX/ 

£ ui(t)  * n hpTi 

i=i  1 i^i  L^i'-x;J 


for  ail  t£n. 


(2) 


Furthermore,  equality  holds  if  and  only  if  there  is  a constant  K 
such  that 


u,(t) 

-T  = K for  all  : 
w . t.  x ) 


This  theorem  is  the  basis  for  condensation  of  posynorm  als . 
Let  t £ n c B^ 


6k 


(t) 


Z u . ( t. ) 

je(k)  J 


c. 

j 


m 

n 

i=i 


j £ \k) 


and 


From  Chapter  t we 
(t)  a "posynomial 


recall  that 

* . ^ " 
term 


is  called  a ‘'posynomtal ' 


ill 


Let  t°  > 0,  t°  £ fl 


and  define 


/ 0 Vt0) 


(5) 


The  theorem  can  now  be  applied 


E 

je(k) 


u (t)  > n 

J je(k) 


(t,t°) 


>) 


>v  . Q 

The  function  g.(t,t  ),  defined  as  the  right  hand  side  of  the 

/ \ 0 

inequality,  is  called  the  condensation  of  g^t)  ^ > the  te™ 

coined  by  Duff in  [2k] . 

The  condensation  of  a posynomial  function  approximates  it  by 
a single  term  posynomial.  The  approximation  has  all  the  properties 
of  a first  order  Taylor's  expansion  of  a convex  function,  as  shown 
by  the  following  lemma. 

Lemma  7.V.  Let  gk(t,t  ) be  as  defined  in  (U).  Then 

(i)  gk(t0>t°)  = gk(t0) 

(ii)  Vgk(t°,t°)  =Vgk(t°) 

(iii )  gk(t,t°)  < gk(t)  for  all  ten. 
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Proof: 


(i) 


Uj'U0)  u^(t°)  *gk(t°) 
Wj(t°)  = 


0, 


0 = gk(t  ) for  all  J € ) 

Uj(t  } 


Thus,  the  condition  for  equality  in  Theorem  7-2  holds. 

Substituting  t^  for  t in  (4)  we  obtain  the  result. 

dgk(t°)  , 0 

(ii)  ~^t—  = “ A,aijUj(t  } 


“i  t i je  <k ) 


dgk(t°,t0) 

ST — 


and  from  ( i ) 

.0  ,0 


_1 

0 


n 


■ / , 0,i 

Uj(M 

0^ 


t"  imrp'’)} 


w,(t  ) 


j£(k) 


ij  J' 


a6.<tU,tU)  , o o 

5+ = ~o  gv(t°)’  E *11*1^°) 

t°  ^ je(k)  J J 


S6k(t°) 


= E aiiui^t0^  = — St 
t°  je(k)  J dt 


(iii)  This  is  just  a restatement  of  (4).  □ 

Figure  1 shows  the  relation  between  g(t)  and  g(t,t^) 


7 -3  Linearization  of  Condensed  Programs 

The  dual  chemical  problem,  stated  in  Cha^t^r  4 is 
Problem  DCP 

m 

Minimize  (-log[gn(t)j  = - £ b ‘log  t } 

i=l 


(1) 
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s(t) 


S“?0I8SS?BS3 

• * • • • 


Figure  1 

Posynomial  Condensation 

g(t)  = l/t  + t + t2 
g(t,l)  = 3-t2/3 


Ilk 


Subject  to  g.(t)  = Z u (t)  <1,  k =1,  2,  ...  , K (2) 
^ j£<k>  J 


and 


te  ijj 


where  u,(t)  is  defined  in  the  preceding  section  and 
J 


Cj  = exp[-Cj],  j = 1,  2,  . . . , n 


Define  the  feasible  set  G of  problem  DCP  by 


G a (t  € rfjlg^t)  < 1,  k = 1,  2,  ...  , K)  (3) 


Let  t £ K+  and  consider  the  totally  condensed  dual  chemical  problem. 


Problem  DCP(t  ) 


m 


Minimize  {-log[gn(t)J  = - Z b ’log  t ) 
u i=l  1 x 


W 


Subject  to  g^tjt*"*)  <1,  k = 1,  2,  , K 


(5) 


t£  ijj 


Define  the  feasible  set  for  the  condensed  problem 


G(t°)  2 (t  £ I^g^tjt0)  < 1,  k - 1,  2,  . ..  , K)  (6) 


Proposition  7.4:  G c G(t°)  for  all  t°  € iP 


Proof:  t £ G implies  gfc(t)  < 1,  k = 1,  2,  . . . , K,  which 

implies  g^tjt0)  <1  by  Lemma  7. 3-  Hence  t £ G(t°).  □ 


A totally  condensed  program  such  as  DCP(tU)  can  be  linearized 


simply  by  taking  the  logarithms  of  all  the  condensed  functions  p . 

&k 


Let 


allt(t°)  * E a v (t°) 
lk  j£(k>  1J  J 


(7) 


ck(t°)  = n 

k J6(k) 


,w^(t0) 


li. 


Vt°) 


(8) 


Here  w(t  ) is  defined  by  (2.3).  Note  that 


Os  . ~ flox  ® aik(t0) 


) = c^t  ) n tj. 


k=l,  2,  ...  , K (9) 


Now,  letting 


z.  = log  t^ 
1 i 


0 , +0 
z,  = log  t. 


and  taking  the  logarithms  of  the  constraint  functions  g (t,t  ) 

K 

in  (2.5),  we  arrive  at 


log  gk(z,z°)  = log  ck(z°)  + Z aik(z°)z;L  <0,  k = 1,  2,  ...  , K 


(10) 
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We  use  the  somewhat  ambiguous  notation  g(t)  = g(z).  It 
will  always  be  clear  from  context  which  function  is  used. 

Note  that  in  the  case  of  a go.  eral  geometric  program,  total 
condensation  also  requires  condensing  the  objective  function.  For 
the  dual  chemical  problem  DCP,  the  objective  function  is  already 
a single  term  posynomia.l  and  need  not  be  condensed. 

With  the  transformations  described  above,  the  problem  becomes 
a linear  program : 


Problem  LP(z^) 


Minimize  -bz 


0.  0. 

Subject  to  A(z  )*z  < c(z  ) 


(ID 


where 


*ik(z0)  ■ aik(z°) 


Ck(z°)  = -log  Ck(2°) 


Notice  that  z is  unrestricted  in  sign 
~ / 0 


Problem  LP(z  ),  approximating  DCP,  has  m variables  and 
K constraints.  For  a typical  chemical  system,  m is  considerably 
smaller  than  the  number  of  species  n.  Typically  m/n  = l/2  or 

even  smaller.  K is  usually  very  small — less  than  10  even  for  large 

0 \ 

problems.  The  matrix  A(z  ) is  thus  much  smaller  than  A,  a 
definite  advantage  in  storage  and  computations. 

Instead  of  working  with  t,  we  could  initially  develop  the 
condensation  with  z = log  t,  that  is,  condense  the  transformed 


i 
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i pfnamp  fpppiip^p  ■ "'^ 


geometric  program  TCP  in  Section  4.4  and  its  corresponding  dual 
chemical  program  DCP.  For  the  transformed  program,  define  the  transformed 
feasible  set 


Z = U € r |gk(z)  < 1,  k = 1,  2,  ...  , K} 


(12) 


Similarly,  for  the  condensed  program  we  have,  based  on  (ll) 


Z(z°)  = (z  £ I!‘  |a(zU)-z  < c(z°)) 


7/  0V  . 0, 


(13) 


Applying  the  strict  monotonicity  of  the  logarithmic  function, 
it  is  easy  to  prove  the  analog  of  Proposition  7.4. 


Proposition  7.5:  Z c Z(zQ)  for  all  z £ 


An  immediate  corollary  is 


Corollary  7.6: 


inf{b'z|z  £ Z(z  ))  < inf[b'z|z  £ Z]  for  all  z £ 


»v  0 

Notice  that  all  the  functions  in  TCP  and  LP(z  ) are  well  defined  and 


■X,  0 


that  Z and  Z(z  ) are  closed,  in  contrast  to  PGP  (Section  4)  and  G. 


Unfortunately  none  of  the  feasible  sets  need  be  bounded.  Even  if  G 


~ 0 ~ 0 

or  Z are  bounded,  G(t  ) or  Z(z  ) may  still  be  unbounded,  causing 


~ / 0, 


possibly  unbounded  solutions  to  LP(z  ).  In  practical  situations  it 
is  always  possible  to  set  upper  and  lower  bounds  on  each  variable, 
e.g., 
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- 

2 < z < z componentwise 


Therefore,  for  practical  purposes  there  is  no  loss  of  generality 
in  assuming  that  Z is  bounded  by  adding  the  bounds  to  the 
constraints . 


7 .4  Convex  Cutting  Plane  Algorithms 

The  reader  may  have  realized  by  now  that  the  idea  of  Dembo's 

method  is  to  generate  a sequence  of  linear  programs  LF  = LP(z  ) 

* 

whose  solutions  will  hopefully  converge  to  a point  z solving  TGP 
or  DCP.  Generation  of  the  sequence  of  problems  is  facilitated  by 
introducing  a "cut"  at  each  iteration,  which  excludes  part  of  the 
preceding  feasible  set  from  the  new  problem.  The  basic  ideas  of 
convex  (or  concave)  cutting  plane  methods  were  formalized  by 
Zangwill  [55 ]■  We  review  here  the  convex  analog  of  his  method. 

The  problem  considered  is 


Minimize  c'x 


Subject  to  g,(x)  <0,  k = 1,  2,  ...  , K. 


The  functions  g (x)  for  all  k are  convex  and  differentiable. 

XV 

There  is  no  loss  of  generality  in  assuming  that  the  objective  func- 
tion is  linear  since  it  can  always  be  replaced  by  a single  variable 
and  added  to  the  constraints . A solution  is  obtained  when  at  some 


iteration  k,  x G G,  where  G is  the  solution  set . 


I 


~ 


Starting  with  a set  IT*-  approximating  0 externally  (that 
is  G 5 we  solve  the  following  problem  LP"'": 

Minimize  cx 
Subject  to  x G U1 

In  general,  at  iteration  k we  solve 


Problem 


LP* 


Minimize  cx 


Subject  to  x 


£ U*. 


A solution  x , to  LP  is  then  tested  to  see  if  x £ G. 

If  so,  the  algorithm  terminates.  Otherwise,  a mapping  A is 

k k k 

used  to  generate  a point  v . Usually  v is  simply  x , i.e., 

k k k 

A is  an  identity  mapping.  Using  v , a half  space  H s H(v'  ) 
is  generated,  having  the  property  that  G c ^ and  x*  $ H*.  The 
set  U*+'*'  is  then  defined  by 


u“+1  . U*  n h* 


k +1  k 

and  we  proceed  to  solve  LP  . The  half  space  H is  generated 
by  a plane  called  a cutting  plane . 

Dembo's  Algorithm  (DA)  replaces  the  objective  function 
g Q(x)  by  a single  variable  xQ  and  the  constraint  gQ(x)  < x . 

In  the  dual  chemical  problem  DCP,  since  the  objective  is  linear  in 
z,  this  substitution  is  not  necessary.  This  fact,  and  other  changes 
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to  his  method  discussed  later,  were  incorporated  in  a modified 
version  of  the  algorithm,  henceforth  referred  to  as  MDA. 


7 • 5 Generation  of  Cutting  Planes 

In  this  section  we  show  how  condensation  techniques  are 
applied  in  DA  and  MDA  in  view  of  the  preceding  section.  Specifically 
we  show  how  the  half  spaces  H^,  defined  by  the  cutting  plane >,  are 
generated.  We  change  our  notation  slightly,  denoting  by  the 

linear  objective  function  of  TGP  and  DCP.  The  linearized  program 
at  iteration  k will  be  denoted  LP  . 

i 

Following  the  notation  of  Chapter  5,  we  denote  the  minimizing 
set  (solution  set)  of  problem  LP  by  D(gQIZ  ),  where  Z is  the 
feasible  set  of  LP  . Note  that  the  superscript  k and  the  subscript 
k are  not  related.  The  former  is  iteration  number;  the  latter, 
constraint  index. 

Lemma  7 . 7 : Let  z £ z"  ))•  Then,  either  there  exists  a k 

such  that  g (z)  >1  or  z £ M(g  |Z)  where  Z is  the  solution  set. 

K 0 

The  proof  is  trivial,  by  definition  of  Z 
<x  k “ X 

Lemma  7.6;  Let  Z = Z(z  ). 

Let  zk  £ D(g0|Zk)  and  suppose  zk  $ D(gQjz). 

Define  L = (k|g^(z)  > l) . 


i 

3 


I 


l 

' 

I 

J 

j 

: 


. 

, 

* 

j 

* 

i 
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For  l £ L let 


Hi(zk)  s (z  € Em|log[gi(z,zk)]  < 0.!  (l) 

Define  Z = Z DH^(z  ). 

Then  the  following  hold: 

(1)  zk+1  = zk 

(ii)  zk?2ktl 
(ill)  Z ~ 2k+1 

Proof i 

(ii)  If  zk  € Zk+1  then  zk  £ H^(zk)  and  thus  g^(zk,zk)  < 1 
would  imply  g^(z  ) < 1 by  Lemma  7-3*  This  is  a contra- 

lr  4-"] 

diction  since  £ £ L.  So  z ^ Z 
(i)  Follows  from  the  definition  of  Z , with  strict  contain- 
ment since  (ii)  holds. 

(iii)  Clearly  Z w Z by  Proposition  7-5-  From  Lemma  7-3  (iii) 

k ~k+  X 

it  follows  that  Z ^ H^(z  )>  hence  ZcZ 
This  completes  the  proof.  □ 

Lemma  7.7  describes  the  solution  test  and  supplies  a stopping 
criterion- -when  g ( z ) < 1 for  all  k.  Lemma  7.8  shows  that  if  a 

Jv  ^ 

stopping  criterion  has  not  been  reached,  a cutting  plane  can  be 
generated  which  meets  the  requirements  of  cutting  plane  methods. 
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Thus,  a new  problem  LP^+^  is  generated  with  a feasible  set  Zk+\ 
Dembo's  algorithm  (DA)  performs  all  its  computation  in  the 
untransformed  variables  t.  In  contrast,  MDA  works  in  the  trans- 
formed variables  z,  for  reasons  discussed  later.  We  show  now  that 
in  the  z-space  the  generation  of  cutting  planes,  and  in  fact,  con- 
densation in  general,  is  actually  a first  order  Taylor's  expansion. 
Recall  that 


g (z)=  E exp[-c  + E a z ] 
k j£(k)  d i ij  i 


We  expand  log[g^(z)]  near  z 


0 


log[gk(z)] 


log[gk(z0)]  + (z-z°)Vgk(z0)-  ~ 

%(z ; 


log{  E exp[-c  + E a .z°A) 

je(k)  J 1J  1 


m 


(z.-z°) 


i i 


i=l  g.  (z°)  j£(k)  1J  " J i 1J  1 


• ( E a.  . exp[-c  i-  E a.  ,z?] ) 


m 


= log  g (z°)  + E (z,-z°)  ( E a w (z0)) 

k i=i  1 1 je(k)  1J  J 


log  gk(z°)  + E Vaik(z°)  - E Vaik(z0)  (2) 

i i 


Here  w.(z^)  and  a,,  (z^)  are  those  defined  by  (2.5)  and  (5-7), 
J ik 

respectively.  From  Lemma  7-5  and  definitions  (3.8)  and  (5*9) 
we  obtain 
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log[gk(z0)]  = log[gk(z°,z0)]  = logCk(z°)  + £ Vaik(zi) 
Combining  with  (2)  we  have 

log  gk(z)  = £ zi'aik(zi)  + loe  ck(z°)  (3) 

The  half  space  H^(z  ) is 

Hg(zk)  = (z  € lf1|Ai  (zk)-z  + log  C^{ zk ) < 0}  (4) 


7.6  The  MPA  Algorithm 

We  assume  that  z is  bounded,  i.e.,  z < z < z componentwise. 
The  bounds  need  not  be  explicit  as  constraints  since  in  the  linear 
programs,  we  can  use  upper  bounding  techniques.  The  lower  bounds 
can  be  eliminated  by  redefining  z and  requiring  nonnega  ..ivity . 

The  problem  solved  is  DCP  described  in  Section  7-3.  We  continue  to 
use  the  term  condensation  although  the  description  is  in  terms  of 
the  transformed  variables  z.  Except  for  the  z-notation  and  the 
fact  that  the  original  objective  function  -bz  is  used,  the  descrip- 
tion below  is  essentially  Dembo's  Algorithm  [20]. 

Algorithm  MPA 

0 ~ / 0\ 

1.  Select  an  arbitrary  point  z and  construct  Z = Z(z  ) by 
total  condensation.  Set  k = 0. 
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2.  Increase  k by  1. 

Solve  LF*  : minimize  -bz 

subject  to  z £ Z 
Let  £ £(  -bz  |Z^) . 

Compute  g (z  ),  l = 1,  2,  ...  , K. 

Let  L = (i  |g/zk)  >1,  i = 1,  2,  ...  , K). 

k i 

If  L is  empty  terminate,  since  z £ £j(-bz|Z). 

3.  Otherwise,  select  i £ L (the  most  violated  constraint  is  chosen). 
Construct  = H^(zk)  by  equation  (5.4). 

~\r  Ir 

k.  Define  Z = Z fllf.  Continue  with  Step  2. 

The  proof  of  convergence  of  this  algorithm  follows  directly 
from  Zangwill's  convergence  proof  [55 1 and  will  not  be  repeated. 
However,  his  proof  depends  crucially  on  the  nesting  property,  i.e., 
the  fact  that 


This  property  can  be  insured,  of  course,  by  retaining  all  previous 
cuts  throughout  the  computation.  With  the  addition  of  cuts,  the  size 
of  the  problem  may  increase  considerably,  and  much  of  the  advantage 
of  the  method  is  lost.  Topkis  [51,52],  and  Eaves  and  Zangwill  [27] 
investigated  conditions  under  which  the  nesting  requirement  can  be 
relaxed.  Topkis  showed  that  if  the  objective  function  is  strictly 
quasiconvex  (ours  is  not),  one  can  relax  nesting  requirements  and 
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Retain - 


k-1  _k 

retain  only  the  active  constraints  of  LP  for  problem  Lr  . 

ing  all  previous  cuts  not  only  increases  the  dimension  of  the  problem, 

it  may  also  cause  ill-conditioning  in  the  LP  matrices,  especially  when 

the  generated  cuts  are  "close"  to  each  other. 

Ihe  nesting  requirement  can  still  be  relaxed  if  one  is  willing 

to  assume,  as  we  are,  that  z is  the  unique  solution  to  LP  . In 

k k+1 

this  case,  maintaining  all  the  active  constraints  of  LP  in  LP 
we  have 


. k+1  k 
-bz  > -bz 


since  zk  $ Zk+1  and  hence  C(-bz|zk+'L)  c:  JQ(-bz|Zk).  In  the  relaxed 

~k+l 

version  of  the  algorithm  only  step  4 is  changed  where  Z is  defined 

by  the  active  constraints  of  Z and  the  half space  H . The  assump- 
k k 

tion  that  each  z be  the  unique  solution  to  LP  is  equivalent  to 
assuming  nondegeneracy  for  each  k. 

Our  experience  showed  that  the  relaxed  version  was  more 
efficient.  The  reduced  overhead  of  maintaining  a large  matrix  more 
than  compensated  for  the  slight  increase  in  the  number  of  Iterations. 

In  the  rest  of  this  chapter  we  discuss  computational  aspects 
of  the  algorithm  and  present  a comparative  study  on  several  test 
problems,  each  solved  by  MDA  and  the  RAND  code  [48]. 
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7.7  Computational  Aspects  of  the  Algorithm 
(a)  Starting  Point 

The  RAND,  DA,  and  MDA  methods  require  starting  points  x^,  t^,  z 

0 0 

respectively.  For  chemical  problems  t and  z may  be  difficult  to 
guess  (although  they  are  not  required  to  be  feasible).  The  overused 
phrase  of  nonlinear  algorithms  "Given  a starting  point  which  is  close 
enough  to  the  solution,  the  method  converges..."  is  true  also  here. 

We  found  that  the  cutting  plane  method  is,  in  general,  more  sensitive 
than  the  RAND  code  to  this  aspect.  The  RAND  code  uses  projection  to 
generate  a feasible  point  x^  if  a nonfeasible  guess  is  given.  If 
no  guess  is  given,  it  generates  a positive  feasible  point  by  defining 

xj ' yj + s 

and  solving  the  linear  program: 

Program  CLP 

Maximize  | 

Subject  to  Ay  +A,|,1  = b 

y > 0 

where  1 is  an  n-vector  of  l's.  This  method,  due  to  Clasen  [12], 
n 

has  several  advantages : 

(i)  It  generates  a positive  vector  x if  one  exists. 

(ii)  It  detects  degeneracies  and  dependent  rows  indicating  infeasi- 
bility when  no  x > 0 exists. 

(iii)  It  does  not  require  starting  guesses. 

(iv)  It  can  generate  a dual  point  t or  z by  using  the  Kuhn- 

Tucker  conditions  and  a least-squares  approach  (see  Section  6.1) 


Despite  these  advantages,  the  initial  point  generated  was  in  many 
cases  very  far  from  optimal,  due  especially  to  incorrect  distribution 
of  species  appearing  in  several  phases.  For  example,  the  bulk  of 
water  is  usually  in  the  liquid  phase,  but  the  procedure  sometimes 
assigns  almost  all  the  water  to  the  vapor  phase  (see  example  below). 

To  remedy  the  situation,  the  following  modified  linear  program  was 
solved: 

Program  LLP 

Minimize  cy  + cl  i 
J n 

Subject  to  Ay  + ln§  = b 
6 > € 

where  c is  the  vector  of  free  energy  coefficients  and  e is  a small 
positive  number  (say  .001)  to  insure  a positive  solution  if  one  exists. 

The  new  procedure  gives  remarkably  improved  starting  points, 
as  shown  by  the  example  below.  The  composition  vector  x and  the 
mole  fractions  x are  given  for  the  two  methods  and  compared  to 
the  equilibrium  values.  Notice  that  the  CLP  solution  has  practically 
all  the  water  in  the  gas  phase  whereas  LLP  correctly  assigns  the  bulk 
of  the  water  to  the  liquid  phase.  Note  also  that  the  dimensionless 
free  energy  F/RT  obtained  by  LLP  is  within  0.2$  of  the  minimum  value. 
Several  problems  were  tested  with  LLP  and  the  example  is  quite  typical 
of  the  improvement  achieved. 
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Example  7.8:  Starting  Points  for  Soda  Pop  Model  (Appendix  A. 4) 


CLP 

LLP 

£ = 0.002 

Equilibrium 

Species 

X 

A 

X 

X 

A 

X 

X 

A 

X 

°2 

5.256 

.039 

5.274 

.056 

5.276 

.053 

co2 

.020 

1.5  do"4) 

6.042 

.064 

6.043 

.060 

N2 

82 . 56O 

.613 

82 . 578 

.879 

82.580 

.826 

h2° 

46.766 

.3^7 

.002 

6.108 

.061 

Gas 

134.602 

93.896 

100.0 

°2 

.020 

.045 

.002 

4.4  do'5) 

co2 

.020 

.045 

.002 

1.3  do"5) 

n2 

.020 

.045 

.002 

3.8  (lo'4) 

H+ 

.051 

.113 

.002 

3.0  do"8) 

oh" 

.020 

.045 

.002 

5.7  (10 '7) 

Cl" 

.080 

.179 

.080 

.0015 

8.0  (10"2) 

.0017 

+ 

Na 

.088 

.197 

.088 

.0017 

8.8  (lO'2) 

.0019 

K 

ao 

-J- 

O 

.108 

.048 

.0009 

4.8  (10'2) 

.0010 

H2CS 


Glucose 


Liquid 


52 . 8c6 
.026 
.002 
.002 
.002 


53  082 


0.995 


46.70 

2.9  (10'2) 
1.8  (10'6) 
4.8  (10’5) 
2.0  (10'2) 


F/RT  -2772.4 


-3122.9 


-3128.9 


' 


.29 


The  (primal)  starting  point  generated  by  LLP  is  used  directly 


to  generate  weights  by  w,  = x..  Unlike  DA,  which  starts  by  condensing 

J J 

at  z , these  weights  are  used  in  MDA  to  condense  and  linearize,  so 
that  is  generated  without  any  point  z^. 


(b)  Bounds  on  Variables 

To  insure  compactness  of  Z , both  DA  and  MDA  require  bounds 

(upper  and  lower)  on  all  variables.  RAND  requires  lower  bounds  on  all 

x.  to  insure  positivity.  Bounds  which  are  too  strict  may  render  the 
J 

problem  infeasible,  whereas  bounds  which  are  too  wide  are  likely  to 
cause  numerical  problems,  and  slow  convergence.  Upper  bounds  on  dual 
variables  can  be  found  using  the  techniques  described  in  Section  6.2. 

In  many  cases  these  can  be  found  by  inspection. 

Finding  lower  bounds  is  a somewhat  more  difficult  problem,  but 
since  the  dual  variables  represent  (theoretical)  concentrations  of  sub- 
species, a bound  which  is  in  the  order  of  x/An  where  A^  is  the 


23 


Avogadro  number  6.02  x 10  , is  certainly  low  enough,  since  A^  is 

the  total  number  of  molecules  in  one  mole  of  any  species. 


(c ) Generating  Cutting  Plar.es  by  Condensation 

A major  advantage  of  MDA  over  DA  is  that  it  works  in  the  trans- 
formed space  (the  z- space)  instead  of  the  usual  dual-chemical  t-space. 
Consequently,  exponentiation  operations  in  DA  are  multiplications  in 
MDA;  similarly,  multiplications  are  reduced  to  additions.  Computations 
in  the  transformed  space  also  eliminate  overflow  and  underflow  problems 
caused  by  large  values  of  the  free  energy  coefficients. 
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(d)  Solving  the  Subproblems 


The  RAND  code  solves  a linear  system  to  determine  a direction 
for  improvement  of  the  composition -vector . DA  and  MDA  solve  LP  problems 
which  are  relatively  small  (in  the  unnested  version)  and  easy  to  solve 
because  the  dual  simplex  method  is  used,  so  that  an  added  cut  is 
actually  an  added  variable.  A feasible  basis  is  available  from  the 
preceding  iteration.  Since  there  are  m dual  variables,  at  most  m 
constraints  are  active  in  each  LP\  The  bounds,  of  course,  are  not 
treated  as  constraints.  The  lower  bound  is  eliminated  by  a trans- 
formation of  variables.  The  upper  bound  is  handled  by  standard  upper 
bounding  methods.  With  the  added  slacks,  the  subproblems  are  of  the 
order  m x 2m.  The  matrix  A is  usually  dense,  so  that  sparse 
matrix  techniques  are  not  applicable. 

(e ) Recondensation 

As  pointed  out  earlier,  an  inherent  problem  in  cutting  plane 
methods  is  ill-conditioning  due  to  accumulation  of  "close"  cutting 
planes.  This  was  apparently  the  reason  for  the  failure  of  the  dual 
methods  in  the  larger  problems  (see  the  computational  results  which 
follow).  Some  ways  to  accelerate  convergence  are  discussed  by 
Dembo  [20] . 

A possible  (untested)  remedy  is  to  "recondense"  the  problem 
after  a fixed  number  of  cuts,  discarding  most  of  the  existing  con- 
straints and  generating  several  new  ones,  by  condensation  of  the 
original  constraints . This  approach  is  analogous  in  some  ways  to 


reinversion  of  matrices  in  LP  codes.-  in  principle,  recondensation 
could  take  place  at  each  iteration,  with  a result  equivalent  to  block 
pivoting.  This  approach  was  tested  and  did  not  prove  fruitful. 


(f ) Selection  of  Cuts 

The  cutting  plane  h\  generated  in  the  k-th  iteration  is  not 
necessarily  unique,  since  any  violated  constraint  g^,  i £ L,  can  be 
condensed.  The  common  practice  of  selecting  the  'most  violated" 
constraint  is  not  necessarily  the  best  policy.  The  situation  is 
somewhat  analogous  to  the  selection  of  entering  column  in  the  simplex 
algorithm,  based  on  the  most  negative  reduced  cost. 

Faster  convergence  (at  the  cost  of  a slight  increase  in  com- 
putation at  each  iteration)  can  be  obtained  by  condensing  the  con- 
straint leading  to  the  largest  change  in  g^  (in  our  case,  an  increase 
in  g^).  In  the  simplex  algorithm  this  can  be  easily  computed,  by 

multiplying  the  ratio  b./a,  . by  the  reduced  cost  c , and  looking 

1 J 

for  the  minimum  over  j such  that  c.  < 0 and  i such  that  a. . > 0. 

J ij 

For  our  purposes,  the  reduced  cost  for  a constraint  can  be  substituted 
by  its  lagrange  multiplier.  Since  the  latter  is  not  known  either,  an 
approximate  multipl:.er--that  of  the  condensed  constraint- -can  be  used. 
This  multiplier  is  available  only  if  the  constraint  was  linearized  in 
a previous  iteration  and  not  discarded.  Then  the  multiplier  A.  of 

A/ 

the  most  recent  linearization  of  g^  can  be  obtained  from  the  solution 

of  (LP)  . We  shall  choose  to  condense  constraint  g in  iteration  k 

s 

when 
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[g  (zk  ) - 1]  •%  = max  [g»(zk  ) - 1]  •%. 

s s ^ L £ 1 

We  have  not  tested  this  scheme  in  our  algorithm. 


7.8  Comparative  Test  Results 

The  tables  in  this  section  present  results  of  test  runs  obtained 
with  MDA  compared  to  results  obtained  with  the  RAND  code  [1+8] . Both 
codes  were  compiled  in  FORTRAN  H,  Optimization  Level  2,  and  run  on  an 
IBM  360/67.  All  computations  and  results  are  in  double  precision. 
Execution  times  in  seconds  of  CPU  include  input  but  not  output  time. 

The  version  of  MDA  used  here  did  not  use  automatic  starting  points 
or  bounds.  Recondensation  and  acceleration  methods  were  not  applied. 

The  following  are  assumed,  unless  stated  otherwise: 

1.  RAND's  starting  point  was  generated  by  program  CLP. 

2.  MDA's  starting  point  was  based  on  a rough  guess,  as  described 
in  Section  6.2. 

3.  The  unnested  version  was  used. 

1+.  Stopping  criteria: 

RAND:  RMS  error  in  mass  balance  < 3 x 10  ^ 

RMS  error  in  mass  action  < 3 x 10  ^ 

MDA:  For  all  k,  g^(z)  < 1 + e,  e = 10"^. 

5.  Accuracy: 

RAND's  results  are  judged  somewhat  more  accurate  based  on 
the  smaller  final  errors  and  the  more  stringent  termination  criteria. 
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Two  sets  of  problems  were  solved,  one  including  chemical 
problems  and  the  other  including  posynomial  geometric  programs.  There 
is  some  indication  that  MDA  did  better  with  posynomial  problems  than  with 
chemical  ones,  due  perhaps  to  the  former's  "dual"  structure. 

All  test  problems  and  their  solutions  appear  in  Appendix  A. 

In  addition  to  the  comparison  of  the  codes,  tests  were  run  to: 

(i)  Compare  nested  vs  unnested  versions. 

(ii)  Test  the  effects  of  starting  points  on  convergence. 

Analysis  of  Test  Results 

Table  1 shows  the  CPU  times  for  chemical  problems.  For  small 
problems  MDA  is  about  20$  slower  than  RAND,  but  for  the  larger  problems 
the  RAND  program  is  considerably  faster  than  MDA,  mostly  due  to  MDA's  slow 
convergence,  attributable  to  the  large  number  of  cuts  and  the  resulting 
ill-conditioning  of  LP  . For  comparison,  Table  2 shows  results  obtained 
for  problems  which  are  geometric  programs  in  nature  (the  last  three 
are  artificial).  It  seems  that  in  this  class  MDA  fares  a little  better. 
The  last  two  problems,  one  with  a dependent  variable  and  the  other 
with  a loose  constraint,  demonstrate  the  advantages  of  dual  methods. 

RAND  was  unable  to  handle  the  i ermer,  while  the  latter,  although 
solved,  caused  several  underflow  warnings,  since  the  values  of  x 
and  x corresponding  to  the  loose  constraint  approached  0,  causing 
numerical  difficulties  in  the  computation  of  F(x). 
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TABLE  1 


TEST  RESULTS --CHEMICAL  PROBLEMS 


Problem 

Dimensions 

CPU( 

sec . ) 

m 

n 

K 

RAND 

MDA 

cuts 

Small  Problem  A.l 

2 

4 

2 

0.25 

0.28 

0 

Soda  Water  A. 2 

5 

13 

2 

0.57 

O.89 

23 

Hydrazine  A. 3 

3 

10 

1 

0.45 

0.50 

12 

Soda- Pop  A. 4 

9 

17 

2 

0.82 

2.06 

61 

Respiratoryl  A. 5 

12 

30 

3 

2.07 

5.80 

1501 

Respiratory2  A. 5 

11 

30 

3 

I.65 

4.49 

110 

Plasma  A. 6 

16 

21 

l 

1.34 

9.16 

1501 

Fetus  A. 7 

19 

51 

7 

5.09 

10.05 

1332 

Notes:  0.  The  problem  names  are  followed  by  their  index  number. 

1.  Failed  to  satisfy  the  convergence  criterion  after  150  cuts, 
terminated  within  0.1$  of  optimum. 

2.  Failed  to  . onverge,  terminated  within  0.2$  of  optimum. 
Dimensions ; 

m = number  of  subspecies 
n = number  of  species 
K = number  of  phases. 
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TABLE  2 


TEST  RESULTS-- POS YNOMIAL  GECMETRIC  PROBLEMS 


Dimensions 

CPU( 

sec . ) 

Problem 

m 

n 

K 

RAND 

MLA 

cuts 

Sea  Power  A. 8 

7 

10 

4 

0.87 

1.14 

35 

Reactor  A. 9 

3 

5 

1 

0.32 

0.41 

11 

Condenser  A. 10 

4 

8 

2 

O.63 

O.58 

13 

Stochastic 
Condenser  A. 11 

9 

13 

6 

1.28 

1.15 

22 

Decomposition  A. 12 

10 

13 

3 

1.48 

2.17 

66 

Dependent 
Variables  A.13 

3 

4 

0 

0.43 

ll 

Loose 

Constraint  A.14 

2 

5 

1 

0.562 

0.4l 

12 

Notes:  1.  RAND  program  cannot  handle  dependent  variables. 

2.  Several  underflow  warnings  were  raised  during  computation. 

Dimensions ; 

m = number  of  variables 
n = number  of  posynomial  terms 

K = number  of  posynomial  constraints,  not  including  the 
objective  function. 
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Effects  of  Nesting  Relaxation 


Table  3 compares  the  run  times  of  MDA  with  and  without  nesting. 
As  expected,  relaxation  of  the  nesting  requirement  increases  the  number 
of  cuts,  but  the  reduced  dimensions  of  the  subproblems  offset  this 
loss  by  reduced  time  per  iteration.  The  results  show  that  the  overall 
effect  is  a reduction  in  computing  time. 

Effects  of  Starting  Point 

Four  problems  were  tested  with  no  starting  point.  In  these 

* 

tests  the  program  automatically  assumed  that  t^  = 1 for  all  i . The 
results  are  shown  in  Table  4--compared  to  the  previous  runs  which  used 
a rough  starting  point  obtained  by  inspection  of  the  matrix  and  the 
free  energy  coefficients.  The  results  show  a marked  improvement 
with  the  better  starting  point.  It  is  expected  that  incorporation 
of  the  technique  discussed  in  7.7(a)  to  generate  a starting  point 
will  significantly  improve  the  run  times. 

Convergence  Rates 

Due  to  the  complexity  of  the  mappings  in  the  algorithm,  no 

attempt  was  made  to  find  a theoretical  rate  of  convergence.  Topkis 

analyzes  convergence  rates  for  cutting  plane  algorithms  [51,52]. 

A simple  test  was  made  to  check  the  "practical"  rate  of  convergence. 

The  objective  function  F(x)  was  computed  after  each  cut  H*  . The 

value  at  iteration  k being  F^.  With  the  initial  estimate  F° 

* 

and  the  final  value  F known,  we  computed  the  "a-cut"  k^,  defined  as 
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TABLE  3 


NESTED  VS  UNNESTED  TEST  RUNS 


Problem 

Nested 

Unnested 

cuts 

v 

cuts 

Soda-Water 

0.72 

12 

0.66 

14 

Soda-Pop 

2.62 

53 

2.06 

61 

Reactor 

0.45 

11 

0.41 

11 

Stochastic  Condenser 

1.25 

20 

1.15 

22 

Decomposition 

2.30 

60 

2.17 

66 

TABLE  4 

EFFECTS  OF  STARTING  POINT  (S.P.) 


Problem 

No  S.P. 

Rough  S 

.P. 

CPU 

cuts 

CPU 

cuts 

Soda-Water 

O.89 

23 

0.66 

14 

Hydrazine 

0.62 

16 

O.50 

12 

Soda-Pop 

5.11 

90 

2.62 

53 

Respiratoryl 

11.76 

116 

8.36 

88 

NOTE:  In  these  two  cases  the  nested  version  was  used. 
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The  numbers  can  serve  as  indicators  of  convergence  rate. 

Table  5 lists  the  res’  -ts  for  a = 0.5>  OC  = 0.95  and 

* 

a = 0.99.  In  this  table  F was  the  final  value  reached,  which 
was  assumed  to  be  the  optimal  value.  In  some  cases  where  the 
starting  point  was  particularly  "bad,"  F®  was  taken  as  the  first  value 
after  some  stabilization  was  obtained. 

It  is  hard  to  draw  any  definite  conclusion  from  these  results, 
but  the  need  for  acceleration  techniques  is  evident  as  some  of  the 
problems  with  the  larger  number  of  cuts  show  slow  convergence  at  the 
tail,  with  the  last  1$  improvement  requiring  as  many  cuts  as  the 
preceding  99$  • 
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TABLE  5 


CONVERGENCE  INDICATORS 


Problem 

-F° 

-F* 

k.  _ 
0.5 

*0.95 

k„  „ 

0.99 

k 

1.0 

Soda-Water 

2241.8 

2253.2 

3 

4 

4 

14 

Hydrazine 

42.822 

47.711 

1 

4 

6 

12 

Soda  Pop1 

3012.4 

3127.6 

2 

13 

29 

53 

Soda  Pop 

3012.4 

3127.6 

2 

13 

29 

61 

2 

Respiratoryl 

1826.8 

1835.1 

14 

45 

75 

150 

Respiratory2 

1809.O 

1835.2 

10 

17 

54 

110 

Plasma2 

830.19 

830.46 

46 

133 

145 

150 

_ , 2 
Fetus 

1866.7 

1870.8 

38 

126 

133 

133 

Sea  Power 

8.42 

126.47 

6 

21 

30 

35 

Reactor 

259.49 

334.26 

1 

4 

9 

11 

Stochastic 

383.89 

883.48 

1 

4 

8 

22 

Decomposition 

1.0 

18.25 

13 

32 

46 

66 

Notes:  1.  Nested  version 

2.  Algorithm  terminated  without  satisfying  convergence 
criteria. 
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CHAPTER  8 


TRANSCENDENTAL  GEOMETRIC  PROGRAMMING 

8.1  Introduction  and  Formulation 

This  chapter  presents  an  extension  to  the  theory  of  geometric 
programming  to  a wider  class  of  functions,  namely,  forms  including 
variables  appearing  as  exponents  (or  in  logarithms).  Although  the 
principal  motivation  is  to  extend  the  applicability  of  primal 
geometric  programs,  our  results  extend  also  to  the  dual,  of  interest 
in  chemical  equilibrium  problems. 

This  chapter  departs  from  chemically  oriented  terminology; 
the  primal  problem  here  will  be  the  primal  geometric  program  and  its 
transcendental  extension,  called  a transcendental  program . The 
equivalent  of  the  chemical  equilibrium  problem  will  be  henceforth 
called  the  "dual  problem."  As  will  be  evident  from  the  following 
section,  many  of  the  useful  properties  of  geometric  programming, 
especially  the  unimodality  of  primal  posynomial  functions,  no  longer 
hold  in  transcendental  programs.  In  this  sense,  transcendental 
programs  have  difficulties  similar  to  those  arising  in  signomial  [39] 
or  complementary  geometric  pregrams  [5] — these  are  the  extensions 
handling  negative  coefficients. 
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Since  neither  primal  nor  dual  transcendental  programs  are 
convex  or  convex- transformable,  Kuhn-Tucker  conditions  are  no  longer 
sufficient  for  optimality.  The  resulting  duality  theory  -will  there- 
fore be  concerned  only  with  stationary  points  of  the  Lagrangean, 
which  may  be  minima,  maxima,  or  neither.  Figure  2 shows  a simple, 
single  variable  transcendental  function  of  the  form  of  interest  in 
this  chapter.  It  shows  that  even  in  simple  cases  multiple  local  optima 
may  exist. 

We  consider  positive  functions  in  positive  variables  t 6 if 
and  0 € P.  The  distinction  between  t and  0 is  mostly  for  con- 
venience. The  functions  are  called  posynentials  (posynomial-exponential), 
and  as  with  posynomial  functions,  we  can  speak  of  posynential  terms.  Thus 
the  form  of  a posynential  function  is 


where 


gk(t,0)  = Z P (t)  Q (0)  R (0) 


p/t)  . c 


m 

n 

i=i 


p 

ii 

1=1 


(1) 

(2) 

(5) 


Ve)  = ^ exp[  Vi] 


(«0 


Here  C.  are  positive  constants  a , b , and  d. 

J ij 

numbers  for  i = 1,  2,  ...  , m,  j = 1,  2,  . . . , n,  i 
(k)  is  a subset  of  consecutive  integers  of  N s (1, 


are  fixed  real 

= 1,  2,  ...  , p . 

2,  ...  , n} . The 
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terms  P,(t)  and  Q . ( 0 ) are  posynomials,  and  so  is  their  product. 

J J 

The  variables  t are  called  the  ordinary  variables , while  the  vari- 
ables & appearing  also  in  the  exponents  are  the  transcendental 
variables . 

We  note  that  if  either  d».  = 0 for  all  i and  j or  b = 0 

■*J  "J 

for  all  i and  j,  then  the  problem  reduces  to  a standard  geometric 
program.  In  the  latter  case  this  reduction  is  achieved  simply  by 
letting  tm+J  - exp^). 

The  Primal  Transcendental  Geometric  Program  is 


Program  PTP 


Minimize 

gQ(M) 

(5) 

Subject  to 

gp, ( lo ® ) < 1>  k = 1,  2,  ...  , k 

(6) 

t > o,  e > o 

(7) 

This  problem  can  be 

characterized  by 

(i)  a positive  n-vector  of  coefficients  C 

(ii)  exponent  matrices  A c R , B £ R , D6  R 

(iii)  a partition  of  the  integer  set  N = (1,  2,  ...  , n)  which 

defines  the  sets  (k),  k - 0,  1,  ...  , K. 

Before  examining  the  necessary  conditions  for  optimality,  we 

give  several  examples  of  problems,  which  can  be  stated  as  transcendental 

programs.  For  convenience  we  define  u , ( t, , 0 ) = 

J 

Since  g^(t.,0)  = E u ^ ( t, 0 ) we  shall  be  interested  in  the  form 

of  the  terms  u.(t,0). 

J 


.(t)  Q (e)  R (e). 

J ' J J 


lU 


m a m b 

1.  Let  u . ( t ) = C.[  II  ti1J]>[  II  {log  ti)  1J]  where  ^ > 1 for  all  i, 


i-1 


i=l 


Define  0.  = log  , i = 1,  2,  ...  , m,  then 


b a . .0 

u.(e)  = c.[n  e iJ)-[n  e lj  i] 

.1  j 1 s 


a.  . 

2.  Let  v (t)  C.  n t,1J  (a  posynomial.  term)  and  suppose  that  we  have 
J i 

a function  of  the  form 

r 

g (t)  = exp{  Z v .( t ) } * v ( t ) 

0 j=l  J 3 

Let  0,=  r1"  , v.(t).  The  function  can  now  be  written 
1 j=l  J 

e 

gQ(t,0)  = e 1-vr(t) 


subject  to 


3.  Let 


Z O^-v  (t)  < 1 
j J 


g (t)  = [ Z V ( t ) ] • [log{v  (t))  ] 

° J=1  j S 


(assuming  vs(t)  > !)•  The  function  becomes; 


subject  to 


gQ(t,0)  = Z 0 v (t) 
j J 


[v  (t)]  1 < 1 

b 


e 0-v  (t)  < 1 
s — 

In  the  three  cases  above,  the  program  is  reduced  to  a transcendental 
program . 
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8.2  Necessary  Conditions  for  Minima 

Proposition  8.1:  Let  problem  PTP  be  superconsistent,  i.e.,  having  an 

interior  feasible  point . 

If  (t,9)  > 0 is  a (local)  minimum  for  problem  PTP,  there 
exist  nonnegative  multipliers  b^,  ...  , 6^)  and  ( ’ * • >\) 

satisfying  the  following  conditions: 

A6  - 0 (1) 


(B  + 0D)5  = 0 


£ 5.  = X, 

je(k>  i * 


xo  ' 1 


Here  A,  B,  and  D are  the  matrices  in  PTP,  & and  A are  the  multi- 
pliers (vectors);  6 is  a pxp  diagonal  matrix  with  = 6 . The 

conditions  above  are; 

( 1 ) the  ordinary  orthogonality  conditions; 

(2)  the  transcendental  orthogonality  conditions ; 


(4  the  normality  condition. 


Proof:  We  apply  the  necessary  conditions  to  the  Lagrangean 


L(t,0,M)s  M0gQ(t,6)  + £ jjjJgjM)  - 1]  (5) 

k=l 


Pk  are  the  Lagrange  multipliers,  p s 1.  Let  (t  ,6  ) be  a 
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minimizing  point  for  PTP.  By  Kuhn-Tucker  conditions  there  exists  a 

ft 

H > 0 such  that 

* * * - *.*■ - (7) 

= 0;  i = i,a,...,P  (8) 

/ 

^[g^(t  ,9  ) - 1]  = 0>  h = 1,2,..., K .(9) 

ft  ft 

Since  t > 0 and  6 > 0 ve  have  by  (9) 


, * * 
L(t  ,9 


) > 0 


So  we  can  replace  conditions  (7)  and  (8)  by  the  logarithmic  conditions 


* * 


8 log[L(t ,9  jU  )]  _ 

a log  ti 


a iog[L(t*,e*jU*)]  _ 
a log  9 1 

. * #.  . *.  # 

Let  u j ( t ,9  ) = Pj(t  ) Qj(0  ) Rj(0  ) 


= 0 


a log  L(t*,e#jU*) li 

* ioe  V " L(t*,e*,/) 


K 


t±  k-0  J€(k)  J 0 


= 0 


(10) 
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b [log  hit  ,0  , |j  ) ] 
d log  e£ 

0* 

i 

, * * * 

Ll,t  ,0  ,|i  ) 

(11) 


K 


K 


— £ Mv(  £ b u (t  ,e  ))+ z Z d U (t  ,e  ) 

a*  ,,_n  k ^(j,  ) J v_n  Kn(c/vSJ!J  J 


= 0 


0,  k=0 


k=0  j£(k) 


where  i = 1,  2,  ...  , m,  i = 1,  2,  ...  , p.  Let 


* * VtV) 

^ 8o(t*,e*)  ' 


j = 1,  2,  ... 


(12) 


Note  that  by  (9),  L(t*,0*,M*)  = gQ(t*,0*).  Interchanging  the  double  sum 
and  substituting  in  (10)  we  obtain  the  ordinary  orthogonality  conditions 
(l).  Substitution  of  6 in  (ll)  and  rearrangement  yields 


£ (Lj  * Vn)5)  * o* 


i • 1,  2, 


» P 


in  matrix  form  these  are  the  transcendental  orthogonality  conditions 
(2).  Condition  (3)  is  a definition  for  A^ . Condition  (A)  is 
satisfied  since  ^ = 1.  Note  that  the  transcendental  orthogonality 
conditions  are  linear  in  the  primal  variables  0 . Observe  also  that 
each  of  these  constraints  has  only  one  transcendental  variable.  □ 


Lemma  8.2:  For  any  feasible  (t,0)  and  any  5(0),  A(0)  satisfying 

the  conditions  of  proposition  (8.1) 


gQ(t,0)  > v(&, A)  = 


C, 


?-D„  6 


(13) 


' ik8 


are  the 


Here  e is  the  base  of  natural  logarithms,  and  D 

l-th  rows  of  B and  D,  respectively.  Since  6 = 6(6)  we  can  write 
v(6,A)  = v(6,X,6). 


Proof:  The  proof  is  based  on  the  geometric  inequality 

(Section  7.2)  and  is  analogous  to  the  proof  of  the  main  lemma  of 
geometric  programming  [25,  p.  114] . Clearly,  for  feasible  (t,6) 
and  A 

K \ 

n [g.  (t,e)]  K < i 

k=l  K 


Therefore 


K \ 

g0(t,e)  > n [gj^Ct,©)]  k 


but 


[«*(*, e)]  k> 


k=0 


n 

j£<k) 


Uj(t,Q) 


6 


i 


\ 


(14) 


(15) 


where  xX  = x X = 1 for  x = 0,  and  u (t,6)  = P (t)  Q (6)  R.(6). 

J J U J 

Thus 


Developing  the  right  hand  side  we  get 
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K \ B 6 

{ n \k)  (n  is  *'  exP(eiDi  5)]} 

k=l  i 


but  by  (2) 

TV  ■ -Bi.6 


so  that  the  last  product  is 

B 

n [6*m  exp(-B^  6 )] 

A second  application  of  (2)  shows  0 = - B.  &/D.  b,  which  after 

substitution  in  the  preceding  equations  leads  to  the  result.  □ 

This  lemma  demonstrates  the  major  difficulty  in  establishing 
a duality  theory  for  transcendental  programs.  The  dual  is  not  a pure 
dual;  its  variables  include  the  primal  transcendental  variables,  and 
therefore  it  cannot  be  maximized  independently  of  the  primal.  The 
appearance  of  e in  the  expression  for  the  dual  function  suggests 
the  following  generalization. 


150 


Corollary  8.3:  Let  Problem  PTP(a)  be  defined  by  replacing  e by 

a > 0 in  PTP,  that  is  to  say,  PTP(e)  = PTP,  and  ,PTP(a)  has  terns 
of  the  form 

d 9 . 

Uj(t,e)  = Pj(t)  Qj(e)  • not 

JL 

Then  all  the  preceding  results  in  this  chapter  hold  for  PTP(a), 

provided  that  D = (d.J  is  replaced  by  D = D log  a. 

*>  J 

Proof: 

dii9i 

a J = exp [ log  a‘d£  9£]  = exp[d  0^]  . □ 


8.3  Duality  in  Transcendental  Programs 

Lemma  8.1:  (Main  lemma  of  transcendental  programming).  Let  (t,0) 

be  feasible  for  PTP, and  let  5(0),  A(0)  satisfy  (2.l)-(2.4). 

Then 

gQ(t,0)  > v(6,A)  (1) 

Moreover,  under  these  conditions 


gQ(t,0)  = v(6,A) 


if  and  only  if 


u.(t,e) 


-< 


E0 


(tyeT 


j e (0) 


Ak(e)>uJ(t,0)  (k),  k = 1,  2, 


(2) 

.,K 
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Proof:  Inequality  (l)  follows  directly  from  Lemma  10.2. 

Suppose  gQ(t,0)  = v(&,A).  Then  in  Lemma  10.2  the  inequalities 
(2.14),  (2.15),  (2.1b)  are  all  equalities.  In  particular,  we  have 
g^(t,9)^  = 1 which  implies  A^  = 0 if  g^(t,6)  < 1 for 
k = 1,  2,  ...  , K.  According  to  the  geometric  inequality  (Section  7.2) 


there  are  constants  q , k = 0,  1, 


, K,  such  that 


= nk  u.(t,e)  , 


J c (k) 


(3) 


Summing  over  j £ (k ) we  arrive  at 


\ = \ 6k(M)  , k = 0,  1,  ...  , K (4) 

so  that  t)q  = l/g0(t,6)  and  hence  (2)  is  satisfied  for  j £ (0). 

If  \ = 0 we  have  by  (3)  6^  = 0,  j £ (k)  and  hence  A^  = 0 so 
that  (2)  is  satisfied  for  j £ (k)  and  q = 0. 

If  TJk  > 0 then  A^  > 0 by  (4)  but  then  g^(t,0) 
implies  gk(t,0)  =1  so  that 
we  see  that  (2)  holds  for  q > 0,  k = 0,  1,  2,  ...  , K. 

j 

Now  suppose  that  (t,0)  is  primal  feasible,  (&(0),  A(0)) 
satisfy  (2. 1-2. 4)  and  that  conditions  (2)  hold.  It  follows  from 
the  geometric  inequality  (Section  7.2),  that  for  each  k > 1 
the  inequality  (2.15)  is  an  equality,  since  if  A,  (0)  = 0 then 

K 

6j(0)  = 0 for  all  j £ (k)  and  then  certainly 


^ = nk  by  (4),  and  again,  using  (3) 
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\ 

g.  (t,0)  k = n 

k j£k 


u^(t,e) 


6. 


6j  J 


= 1 


If  A^(0)  is  positive,  then  so  is  6^(0)  ^or  a-*--L  j £ (k)  , k > 1, 
But  then  (2)  implies  that 


g Jt,e)  = £ . u,(t,e)  = — Z 6,  = l 


,el>  J 


{-3 


and  consequently 


i = n 

j£(k) 


u,^( t,8 ) 


lJ 


'j  J 


Thus  the  right  hand  side  of  (2.15)  is  always  1 for  k > 1.  For 


gQ(t>0)  the  geometric  inequality  leads  to 


gn(t,e)  = n 

0 j£(o> 


Combining  the  observations  above,  we  get 


6j 


>6, 


g0(t^)  = 


n /u(t,e)\&J] 

n I 


J-H  uj 


K \ 

n v K 

k=0  * 


(5) 


and  since  5(0),  A(0)  satisfy  conditions  (2. 1-2. 4),  the  right  hand 
side  equals  v(6,A).  □ 
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This  lemma,  which  is  analogous  to  the  main  lemma  of  geometric 

programming,  does  not  provide  an  analogous  dual  program.  Note  that 

maximizing  v(6,A)  subject  to  (2.l)-(2.4),  5 > 0,  A > 0,  8 > 0, 

does  not  provide  any  solution  to  the  primal  since  8 may  be  primal- 

infeasible.  Moreover,  even  if  8 is  primal  feasible,  the  solution 

may  be  a local  minimum,  a stationary  point,  or  even  a local  maximum. 

Under  more  restrictive  hypotheses  we  can  achieve  a more  tangible  result. 

First  we  wish  to  verify  that  for  any  fixed  8 for  which 

both  primal  and  dual  problems  are  feasible,  the  duality  relations 

reduce  to  the  usual  geometric  programming  duality.  To  see  this,  let 

8 be  fixed  so  that  for  every  1 < j < n,  Q.(0),  R,(0)  is  fixed. 

0 J 

Define  C.  = C Q.(0)  R (0)  and  note  that  PTP  is  now  a geometric  program, 
J J J J 

PGP  (Section  4.1),  with  a new  coefficient  vector  C.  Looking  at  the 
dual  function  v(5,A): 


v(6,A)  = 


n 

n 

j=l 


u (t,e)\  J 

„U...  - - 
6 . 

J 


' n A' 

k=l  * 


I 

fi[>' 

5J- 

- 

II 

J 

. 

n (Qi(e)  R«(e)) 

j=i  J J 


j 


K A' 


n 

j 


c.iQj(e)  Rj(0)\  j 


'j 


k K 


(6) 


o,\6j 


n 


L M J J 


f A 


k 
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Since  0 is  fixed,  dL/50  = 0,  and  the  optimality  conditions  on  9 
in  the  dual  are  now  redundant.  The  remaining  constraints  are 
the  usual  dual  geometric  program  constraints,  so  that  the  dual  problem 
is  the  geometric  programming  dual,  with  the  new  coefficient  vector  C. 

Our  next  theorem  is  the  analog  of  the  duality  theorem  of 
geometric  programming.  To  state  and  prove  it,  we  first  state  the  dual 
transcendental  program . 

Program  DTP 


n ci\  •> 

min  max  v(B,A ,0)  = ii 
0 A L J=1  \ *3/ 

’ K V r bii 

n \ n in  l) 

>1  k J [}  t 1 

5 

J 

subject  to:  0 is  primal  feasible 

(7) 

AB  = 0 

(8) 

(B  + 0D)&  - 

0 

(9) 

Z 6 = < 

JS(k)  J 

f A^j  k = 1,  2,  ...  , K 

Jl  , k = o 

(10) 

6 > 0, 

0 > 0 

(11) 

Theorem  8.5;  Suppose  the  primal  program  PTP  is  superconsistent  and 

* * 

attains  its  constrained  minimum  at  a feasible  point  (t  ,9  ).  Then 
(i)  The  corresponding  dual  program  DTP  is  consistent  and  has  a 
feasible  solution  (&  ,A  , 6 ). 

(ii)  gQ(t  ,0  ) = v(B  ,A  ,0  ). 

^iii)  the  relations  ().2)  hold  at  (t  ,0  ,&  ,A  ). 
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Proof; 


(i)  By  Proposition  8.1  the  dual  is  consistent.  We  now  show  that 
* 

the  vector  & defined  by 


t . * 


Uj(t  ,O/g0(t  ,6  ) 


* U.(t  >d  ) 

1 “k^vj 


<0) 


i e 00 


is  dual  feasible  where  ^ are  the  Lagrange  multipliers  whose  existence 

# * 

is  assured  by  Proposition  8.1.  Clearly  5j  > 0 and  by  hypothesis  9 > 0. 
Thus 


j iJ  J gQ(t*,0*)  k=0  ^ j€*i)  iJ  J 


?.a,4uf(tV), 


i = 1,2, . . . ,m 


By  Proposition  8.1  the  right-hand  side  equals 


o log[L(t*,0  ,u*)] 
8 log  t± 


# # # 

which  must  be  0 when  (t  ,0  , ^ ) are  the  optimal  values.  Similarly 

J + 


K 


(A  A .A  V.,(t  )] 


. (t°,»°)  1=0  J£(k)  -«'J' 


K 


+ 9]  z £ d u (t*,0*)]}=  ° by  (2.11) 

k=0  ,1€  (k)  J 


J€(k) 
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It  remains  to  define 


Z 5 
j€U>  J 


/ * * \ 
g0(t  >d  ) 


to  complete  this  part  of  the  proof. 

* * * 

By  Lemma  8 4 the  relations  between  5.  and  u (t  ,6  ) imply 

J J 


, * *,  . * * *, 
gnU  >d  ) = v(6  ,A  ,0  ) 


It  follows  from  Proposition  8.2  that 

v(&  ,A  .0  ) = max  v(&,A, 6 ) (12) 

&,  A 

subject  to  (8)-(ll)  with 

* 

© = © . 

Now  suppose  that  there  is  another  pair  (t^,  0#)  satisfying  the 
necessary  conditions  for  a minimum  but  with  ® )• 

Using  arguments  similar  to  those  above,  we  infer  the  existence  of 
6^,  A^,  satisfying  the  d .1  constraint-  and  conditions  (4.2)  so  that 
we  must  have 

g0(t#,e#)  - v(b#,A#,e#) 

= max  v(B ,A,0#) 

5, A 

subject  to  (8)- (.11)  with 

0 = 0^ 
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Thus  max  v(5,A,0")  > max  v(8,A,0  ).  Furthermore,  given  0 , we 

. * 

have  shown  earlier  that  min^  ) subject  to  t > 0 and 

* 

g (t,0  ) < 1,  is  an  ordinary  geometric  program  having  a corresponding 

K 

dual  problem 

max  v(B,A,0  ) 

5,  A 

subject  only  to  (8),  (10),  (’ll)  without  the  constraint 

(B  + 0D)5  = 0 

. . * 

It  follows  from  (12)  that  if  0 is  optimal  then 

max(v(5,A,e*)|(8)-(ll)}  = max(v(5,A,0  ) |(8),(l0),(ll)} 

6, A 6, A 

Hence  .(9)  can  be  eliminated,  as  long  as  0 is  primal  feasible. 
Therefore^  one  can  view  the  primal  problem  as 

Program  PTP 

minlmin  g(t,0)}  '(13) 

0 t 

gk(t,0)  < 1 (14) 

t,0  > 0 (15) 

In  this  representation  the  dual  is 
Problem  DTP1 

min  max  v(b,A,0)  (l6) 

0 6, A 
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subject  to 


AS  = 0 


(17) 


E - 1 (18) 

j£<0>  J 

E 5 = \ , k = 1,  2,  ...  , K (19) 

je<k)  J * 

8,X  > 0 (20) 

Q primal  feasible  (21 ) 

Since  we  have  shown  redundancy  of  (9)>  DTP1  is  equivalent  to  DTP 
when  9 is  primal  feasible.  This  completes  the  proof.  □ 

The  dual  problem  DTP  is  unsatisfactory  in  its  current  form. 

The  reader  may  justly  ask  how  the  condition  of  primal  feasibility 

could  be  incorporated  in  a dual  problem,  without  attaching  the  whole 

primal  constraint  set  to  the  problem.  The  restriction  on  9 serves 

* * 

to  insure  global  minimum  of  g(t  ,9  ),  but  the  theorem  sacrifices 
much  in  elegance  to  obtain  it.  In  some  situations,  9 may  not  be 
severely  restricted  by  the  primal  constraints,  for  instance,  in 
primal  unconstrained  problems.  In  this  case  conditions  (7)  can  be 
ignored. 

We  formulated  the  dual  problem  with  the  variables  9 appear- 
ing in  the  objective  function.  It-  is  easy  to  express  9 explicitly 
in  terms  of  B and  D by  (9)  and  then  substitute  in  v to  obtain 

an  apparently  "pure"  dual.  Specifically,  we  find  9 = -B  &/D  & 

X X • X • 
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as  shown  in  the  proof  of  Lemma  8,2,  and  we  can  use  Equation  (2.13). 
The  dual  problem  obtained  this  way  is 


Program  DTP2 
max  v(b,X) 


- 

fcj 

v 

r \ "I 

- 

! B.  b ' 

B,.6' 

n 

In  X k 1 

n 

£■ 

b , 

1 , k j 

eD . & 

.J 

il 

Lk  J 

& 

i i 

subject  to  Ab  = 0 


£ 6 = X,  , k = 1,  2,  , K 


£(k> 


J « 


V1 


B > 0 


(22) 

(23) 

(24) 

(25) 

(26) 


This  problem 
be  undefined 
in  (22).  If 
define 


is  not  equivalent  to  DTP.  In  particular  v(b,X)  may 


even  for  positive  vector  5 due  to  the  last  product 

* * 

a finite  solution  6 , X exists  to  DTP2,  one  can 


* 

d 


B, 

it 


*• 

5 


* 

0 


All  conditions  fcr  a stationary  poin’  ir.  the  lagrangean  L(t,0) 
exist,  except  (perhaps)  feasibility  of  0.  At  any  rate,  these 
are  necessary  but  not  sufficient  for  a minimum  of  the  primal. 


l6c 


8.4  Some  Properties  of  Transcendental  Programs 


( a ) Recovering  primal  variables  from  the  duaL  solution 

As  noted  earlier,  0 need  no*  appear  explicitly  in  the  dual, 

although  its  implicit  value  appears  in  the  dual  objective  function 

* * 

and  must  be  primal  feasible.  If  we  have  a solution  6 , X to  the 

* 

dual,  0 can  be  easily  recovered  by 


* bJ  6 
*1  " * 

*•  • 


Once  0 is  known  it  can  be  incorporated  into  the  constants  C,,  and 

J 

the  problem  reduces  to  the  ordinary  geometric  programming  primal-dual 
relations,  discussed  in  Section  4.3. 


(b)  Degrees  of  difficulty  in  transcendental  programs 

Although  the  constraint  set  of  DTP  is  more  restrictive  than 
that  of  the  ordinary  dual  geometric  program,  the  ordinary  orthogonality 
constraints  Afc  = 0,  ar.d  the  normality  constraints  which  relate  to  the 
nontranscendentai.  variables,  s*ili  determine  the  dimension  of  the  dual 
space.  Indeed  if  r.  =-  m+1  and  rank  A = m,  these  equations  will  have 
a unique  solution  (which  may  or  may  not  be  feasible  even  if  5 > 0). 

We  define,  therefore,  the  degree  of  difficulty  ir.  the  same  way  as 
fcr  ordinary  geometric  programs,  dd  = n-m*l.  If  dd  = 0,  and  the 
dual  is  feasible,  the  solution  is  unique.  This  striking  feature  is 
demonstrated  by  the  following  example. 


l6l 


; 

j 

i 


Example  8.6: 


PTp 


min  gQ(t,0) 

= 2-t1’ti+  + t^1-t2-01-02,exp(e1)  + 1 ^•t5-ti+1-01>exp(-3ei+e2; 


subject  to 


gl(t,e)  = t^ye'  + t2-t^ 

•02  -exp(-202)  < 1 

t,  0 > 0 

The 

dual  constraints  are; 

1. 

&1 

+ 62  + 

= 

2. 

V 

61 

&2 

O n 
II 

3. 

V 

62  - 4&3  t 

\ * V 

4. 

t : 

6 

& = 

3 

3 

5 

5- 

V 

51 

63 

= 

6. 

V 

(1  + ex)b2  + (1  - 3e1)&5  - 

&4 

7. 

V 

- (1  + 202)&5  = 

1 normality 


ordinary 
) orthogonality 


0 | transcen- 
dental 

orthogonality 


The  first  five  equations  have  the  unique  solution 


* (l  l l 1 !\ 

0 _ '>4’  2’  k’  4’  4' 


Substituting  into  equation  (6),  we  have 


= 2 
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From  equation  (7)  we  have 


The  dual  objective  function  is 


I 


Hence 


t,  - 8V8  e-n/8 

tx  . a'1'8  e11/8 

H .is1/1*  e-3.1*/(8-V8eU;8|  , i8V8e-17.8 


•l  -j  : 

? _ f,  e ~ ■— 

2 y 2 


, -2  „ q1/8  -27/8 
= 2e  tg  = 2*8  ' e 


Thus 


t*  , (8-1/8eU/8>8l/85-ll/8;2.8l/8e.27/8j 


To  verify  that  this  point  is  not  a maximum,  let  t^  increase  without 
limit . 


8 t Condensing  Transcendental  Programs 

In  much  the  sane  way  tr.at  condensation  was  applied  to  geometric 
programs  in  Chapter  9>  it  can  te  applied  also  to  transcendental 
programs.  In  fact,  the  proof  of  Lemma  8.2  already  showed  a condensation. 
In  general,  let 

>d  = Z u (t,e 
J€(k>  J 

and  let  (t C,6°)  be  some  positive  point.  Define  the  weights 
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(1) 


o_  ,t0fl0,  _VtV) 

w.  = w (t  ,e  ) = — ___ 

J J gk(tu,eu) 


U (t,0)\  J 


The  right-hand  side  is  a single  term  pcsynential  having  new  exponents 


a,.  ” a (t°,e°)  = Z a 

lk  ik  J4{k)  U j 


*tk  3 fVt°’e°>  ’ ^ 


r/k  * rJk(t°'e0) 


and  a new  coefficient 


v 5-“Vl  ■ |f 


g k<  t,e,t°,e°)  - ck  n t ik  n n exp(e  Tik) 

i i i 


Taking  logarithms,  we  obtain 


ek  , log  ck  * Z oik  log  * Z pfk  log  el  * Z rfkaf  (7) 

1 H l 


This  form  is  a linear- Logarithmic  function. 


A condensed  transcendental  program  PTP  can  thus  be  viewed  as 


a General  Linear-Logarithmic  problem.  In  this  form,  however,  we  do 
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not  have  the  product  Q In  9 which  characterizes  Clasen's  Linear 


Logarithmic  problems  [12].  It  can  be  easily  verified  that: 

(i)  ik(t°,e0,t°,e0)  = gk(..°,e0)  (8) 

(ii)  Vg  {t°,e0,t°,e0)  = vgk(t°,e°)  (9) 

k K 

(iii ) gk(t,,0,tC,0°)  < gk(t,e)  for  all  (t,0)>0  (10) 

The  proof  is  analogous  to  that  of  Lemma  7*3- 


8.6  Solving  Transcendental  Programs  by  Condensation 

In  principle,  we  could  apply  the  same  ideas  used  in  the  cutting 
plane  method  of  Chapter  9 to  transcendental  problems.  Regrettably,  our 
subproblems  are  not  linear,  but  are  linear- logarithmic  (or  linear- 
exponential).  Trying  to  'linearize'1  (7)  by  defining 


z . = In  v 

i l 


ye  s ir.  et 


leads  to  the  function 


y 

5 z Vi  " ? + ? Ik  e + 1o« 

x /.I  a 


Noting  that  e is  a convex  function  on  1,  we  can  still  linearize 
by  observing  that 
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> / * (y  - y0)ey° 


Letting  y^  = log  0 ^ the  condensed  linearized  form  of  the  constraint 

g (t,6)  < 1 becomes 
K 


r"(^y,t°,e°)  = Z alkzi  * Z ( e<k  * rtk)vt  + \<° 

* jL 


where 


Y = 6°  Y 
rek  £ r£k  ' 

5k“ ln  \ ^rikei(1 


log  9°z) 


If  x > 0 for  all  i then 

JtrL 


g(Zyy)  > g(z>y,t°>9°)  > gL(z,yyt°,0O) 


and  all  the  cutting  plane  arguments  remain  valid,  so  this  could  serve 

as  a solution  method.  If  y ^ are  not  nonnegative  (and  recall  tha-' 

these  are  functions  of  the  condensation  point  and  not  constants!), 

then  the  right  inequality  above  need  not  be  valid,  in  fact,  we  may 
~L  v. 

have  g g.  The  algorithm  in  this  case  no  longer  guarantees  the 
nesting  property  and  the  "cuts"  need  not  exclude  the  previous  point, 
so  that  the  metnod  is  no  longer  applicable.  If  d„  > 0 for  all  11 
and  j then,  by  using  the  transformation  z.  = log  t.  and  y = log  6 

1 1 Ju  At 

we  obtain  a convex  transformed  transcendental  program  with  transformed 
posynential  terms  of  the  form: 
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Uj(t,y)  = C.  exp[E  a^Zj)  exP^£  b^y^}  exp(£  d^e  *) 


-c3  expl?  Vj  +^[Vj  + V 11 


In  this  case  condensation  at  any  (t,0)  > 0 will  yield  nonnegative 

y , so  that  linearization  is  possible. 
i/K 

It  seems  that  in  general,  the  primal  transcendental  program 
is  easier  to  solve,  as  the  transformed  problem  is  well  defined  for 
all  y,  z.  The  dual  problem  may  be  ill-defined  even  for  positive  6 
The  resulting  computational  difficulties  may  offset  any  advantages 
offered  by  linearity  of  the  dual  constraints. 


8.7  Nonideal  Systems  and  Transcendental  Programs 

Most  of  the  computational  and  theoretical  treatment  of  chemical 
equilibrium  presented  so  far,  was  confined  to  ideal  system.  Reviewing 
some  notions  of  Chapter  2,  recall  that  the  (dimensionless)  free  energy 
for  a general  system  is 

F(T',P,x)  = £ x m (T,F,x)  (1) 

J ° J 

where  jji.(T,P,x)  is  the  chemical  potential  or  the  partial  molar  free 
J 

energy.  For  fixed  temperature  T and  pressure  P we  obtain  from 
(2X5)  after  dividing  by  RT 

Hj  = + lo6  a (2) 
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In  principle,  all  of  the  nonideality  in  the  system  is  expressed  in 
a . the  activity  of  species  S..  In  the  ideal  case  a.  - x..  In  the 

y ^ j j ,i 

nonideal  case,  without  loss  of  generality,  the  nonideality  can  be 

lumped  in  the  activity  coefficient  y .(x)  by  defining 

J 

a = T . ( x ) • x (.'/) 

J J J 


By  this  definition  of  y (x),  the  free  energy  of  the  system  becomes 


F(x)  = Z x.(c.  + log  y.(x)  + log  x ) 
j J J J J 

J 


(W 


The  term  log  y.(x)  represents  the  deviation  from  ideality. 
J 

It  would  be  nice,  theoretically,  if  F(x)  retained  its  cherished 

properties--homogeneity  of  degree  one  and  convexity.  In  practical 

applications  this  need  not  be  the  case,  as  y\(x)  is  sometimes 

expressed  by  empirical  relations.  Still-,  in  the  majority  of  cases 

at  least  the  homogeneity  is  preserved.  Recall  that  y.(x)  must 

J 

be  homogeneous  of  degree  0 to  maintain  homogeneity  of  F,  since  it 


is  multiplied  by  x . • This  suggests,  as  is  indeed  t-rue  by  chemical 

arguments,  that  y,(x)  is  a function  of  concentration,  independent  of 
J 

the  total  amount  of  one  specier  or  another.  That  is  to  say, 


r (x ) = r,(S) 

J J 


Several  relations  of  this  type  appear  in  the  chemical  literature;  some 
of  the  useful  relations  are  summarized  by  Prausnitz  et  al.  [42] . 
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I 


For  a vapor  phase,  perhaps  the  most  commonly  used  correction 


is  by  the  virial  equation.  With  it,  the  correction  takes  the  form 


108  rj(x)  ■ ; pipji + ^ * • ' ' <5) 

Here  v is  the  molar  volume  of  the  phase,  (3,.  and  y are 

Jl  J1K. 

empirical  constants,  usually  determined  by  complex  relations  to  other 
parameters . 

For  liquid  phase,  the  Wilson  equation  is 


log  rj(x) 


1 - log[£  Xj.Ajj.] 


Z 

i „ Ai 


where  again,  the  constants  are  determined  by  complex  relations. 

We  shall  see  that  the  dual  transcendental  program  leads  to 

corrections  of  the  form  y.(x)  - y.(x).  The  exact  form  of  y (x) 

J J J 

and  its  meaning  depend  on  appropriate  choice  of  coefficients  in  the 
matrices  B and  D of  the  transcendental  program.  In  this  sense, 
transcendental  programs  are  a generalization  of  the  ideal  chemical 
equilibrium  problem.  We  consider  here  program  DTK?  of  Section  8.5. 
Taking  the  logarithm  of  the  dual  function,  replacing  6 by  x and 
A by  x to  conform  with  chemical  terminology,  we  find 
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F(x)  = - log  v(x,x) 


= E X.[-log  C.  + log  X ] - Z X log  X 
. J J J k K K 


, J 


J 

1 V 1 

n 

z 

Z (B  x)  log 

j 

£ 

1 '■ 

iM 


- log  C.  + log  x,  - Z b log 


V 

Di.x 


CJ  + log  *j + f bo loe  \ - V1 


Here  c.  = -log  C +Zb.  B , D are  the  £~th  rows  of  B and  D, 
J J SC  SC  j • 

respectively.  With  c.  taken  as  the  new  free  energy  coefficient,  the 

J 

sum 


I.  d x. 

z b/i log I"  'z~hAi 

l zj  Vj, 


(7) 


is  the  logarithm  of  the  'correction"  term  y .(x),  accounting  for  non- 

J 


ideality.  Thus 


r4(x)  - n 


£.  d x 

_si Li- 1 


£.  b x . 
J i.J  J 


(8) 


The  bracketed  term  must  be  positive  for  y.  ( x ) to  be  well  defined. 

O 


i n 


Based  on  the  structure  of  B and  D we  can  now  examine 
several  cases.  Note  that  B and  D have  n columns  and  p rows 
where  p can  be  arbitrarily  chosen  (equivalently,  the  primal  program 
has  p transcendental  variables).  We  shall  choose  p to  yield 
meaningful  expressions  for  ifj(x). 


1.  p = K (the  number  of  phases) 
Let 

b, 


ij 


For  example 


1 J c (J> 

0 otherwise 


110 0 

0 0 1110.  . . 0 
0 . . . 0 1 0 . . 0 
0 . . . .01111 


Each  row  represents  a phase.  There  are  l's  in  the  row  in 
corresponding  to  species  in  the  phase.  Then  for  species 
we  obtain  by  expanding  (8) 

K. 


rj(x) 


■Lh  tA 

^it(k)  Xi 


z 

i€(k> 


d £. 
ki  l 


the  locations 

S in  phase 
J 


1 


i 

i 


This  introduces  for  each  phase,  a single  correction  factor  y for 
all  species  in  the  phase.  The  factor  depends  on  the  concentrations 
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of  all  species  in  that  phase.  We  have  been  unable  to  find  such  ar. 
example  in  the  chemical  literature  although  one  can  conceive  of  a 
situation  where  this  will  be  approximately  the  case,  e.g.,  when  a 
polar  solvent  "dominates"  the  activities  of  all  species  in  a 
solution . 


2. 


P = n 
Let 


(number  of  species) 

/ 


b 


B is  an  n x n diagonal  matrix 


l = j 

otherwise 


Let 


d-«j  = 


k(i)  = k(j) 


0 otherwise 

D is  an  n X n matrix  with  block  angular  form.  Each  block 
corresponds  to  one  phase.  For  example 


'dl  dl  •••  dl 
d2  d2  ' ' ' d2 

d d . . . d 

n,  n n. 

i l.j 

0 

o 

0 

d , . . . d , . 

n^+1  Hj+1 

o 

Q.  • • t Cl 

n2  n2 

o 

0 

• • • a.  +1 
“2  * “2  " 

d d 

- 

n3  ' n3 
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f 


Let  S £ 4>  and  expand  (8)  with  the  given  B ana  D to  get 
j k 


r .(x)  = — 

J d 

L J 


-b  x. 


* w 3 


where 


P.  = - b ,/d. 
J jJ  J 


n , = -b  . . 
'j  JJ 


By  appropriate  selection  of  dy  P and  r\  are  independent.  Vj(x)  is 
then  a two  parameter  correction.  However,  taking  the  logarithm 


log  r.(x)  = log  p + n.  log  x 
J J ■**  J 


we  see  that  log  p . can  be  lumped  with  c . and  the  remaining  part 
J J 

is  a linear  term  in  log  x..  This  is  perhaps  the  simplest  and  most 

J 

useful  correction  as  it  corrects  the  "ideal"  log  x,  to  (rj  +l)*log  x 

J J J 


i.e., 


F{ x ) = E x . . c . + n . log  x . ) 
' j J J 'j  J 


where  c includes  log  p.  and  t)  . = rj . + 1.  It  is  relatively 
J J J J 

easy  to  construct  empirical  corrections  of  this  type.  In  some  sense, 
this  is  a "first  order"  correction  for  nonideality. 

In  a similar  fashion,  by  selecting  more  complex  structures 


for  B and  D,  more  elaborate  corrections  can  be  obtained.  Notice 
that  in  our  corrections,  y.(x)  is  homogeneous  of  degree  zero,  as 

' d 


required  by  theoretical  considerations.  We  do  not  propose  here  a 
solution  technique,  nor  do  we  imply  that  transcendental  duality  promises 
a new  theory  of  chemical  duality  for  nonideal  system, analogous  to 
the  theory  for  ideal  systems.  Still,  the  examples  above,  and  potential 
applications  of  transcendental  programs  in  engineering  design,  suggest 
that  the  first  steps  developed  here  merit  further  study. 


APPENDIX  A 


TEST  PROBLEMS  AND  SOLUTIONS 


Seven  chemical  problems  and  seven  geometric  programming  problems 
are  presented  here  for  reference.  Most  of  the  problems  were  found  in 
the  literature  but  in  some  cases  these  were  slightly  modified.  All  the 
solutions,  except  as  noted,  are  those  obtained  by  RAND's  Chemical 
Composition  Code  [48].  MDA  solutions  vary  slightly  and  are  somewhat 
less  accurate. 

Data  contained  for  chemical  problems: 

1.  Name  and  source  of  problem. 

2.  Dimensions:  m = number  of  subspecies, 

n - number  of  species, 
k -■  number  of  phases. 

* . 

3.  F /RT  = minimal  dimensionless  free  energy. 

4.  Subspecies  table: 

i = serial  number, 

= name  {formula)  of  subspecies, 
b^  = right  hand  side  (moles  of  B^), 
z^  = optimal  multiplier  (dual  variable). 
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5.  Species  table: 


k 

j 

S . 
J 

c . 
J 

formula 

*■ 


phase  number, 

species  number  (sequential), 

species  name  (formula), 

free  energy  coefficient, 

composition  of  S.  by  subspecies, 
J 

equilibrium  composition  (moles). 


After  each  phase  its  name  and  its  total  number  of  moles  ar.  ]ist.e 

Problems  A.8-A.14  are  geometric  programming  problems.  Data 
for  these  problems  includes; 

1.  Name  and  source  of  problem. 

2.  Dimensions:  m = number  of  variables, 

n = number  of  terms, 

k = number  of  constraints  (without  objective  function 


3.  Terms  table; 

k = constraint,  number  (0  = objective  function), 

j = term  number  (sequential) 

C.  = coefficient 
J 

Exponents  - variables  and  their  exponents  in  term  j, 

* 

a.  = optimal  value  of  the  term. 


Each  constraint  is  followed  by  the  value  of  its  multiplier  at  the 

optimal  solution. 
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4.  Variables  table: 


i serial  number, 

symbol  ~ symbol  of  variable  i, 

Low  bd.  = lower  bound  used  in  MDA  computations 

Upper  bd.  = upper  bound  used  in  MDA  computations 
* 

t.  = optima,  value. 

i 


Note:  idle  notation  2.5(iO  ^ ) means  2.5  x 10 


A . 1 SMALL  PROBLEM  (Clasen  [14]) 


I 


t 

! 


i 

( 


Dimensions;  m = 2, 


n = 4,  k = 2 


F*/RT  = -5-6755 


1 

i 

B. 

b. 

* 

Z . i 

1 

1 

1 

1 

HI 

3-5 

-1.1032 

2 

R2 

4-5 

-0.4052 

k 

— 

J 

s. 

J 

Formula 

4 

x . 

.] 

1 

1 

Cl 

0.7 

R2 

0.8540 

2 

C2 

-0.7 

HI 

1.6795 

PHASE1  x = 

2.5156 

2 

3 

C3 

0.0 

R2 

3 . (iCbCi 

k 

C4 

0.0 

R1 

1 .8205 

PHASE2  xr  -- 

2 

5 . 48c>5 
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A. 2 SODA-WATER  (De  Haven  [19]) 

Dimensions:  m - 5>  n = 13,  k = 2. 

f*/RT  = -2253.21 


0.6469 

0.2602 

3-7058 

55-8103 

55.8103 


-IO.85 
- 7.69 
■11.52 
■36.61 


formula 


GAS  PHASE 


-11.8027 

-24.8080 

-14.5820 


0.6468 

O.2588 

3.7054 

0.3050 


4.9159 


5 

°2 

0.0 

°2 

1.42  (ic’4) 

6 

cvt 

0 

0 

.0 

C°2 

1.34  /10’3) 

7 

N2 

0.0 

N2 

4.15  do’4) 

8 

H+ 

0.0 

H+ 

2.58  (10’3) 

9 

oh" 

0.0 

oh’ 

9.34  (io‘io) 

■+* 

10 

h2o 

-39.39 

H + OH 

55.5053 

11 

HCO’ 

-20.86 

co2  + oh" 

2.58  do’3) 

12 

H2C°3 

-33-61 

C02  + H+  + OH" 

4.13  (I0"b) 

13 

co" 

6.73 

C02  - H+  + OH’ 

5.78  (io'n) 

LIQUID  PHASE 


55.5073 


Dimensions;  m = 3,  n = 10,  k = 1 
F*/RT  = -V7.7611 


l Dantzig  [54]) 

; = 1 

b. 

1 

* 

z. 

1 

2.0 

- 9.7851 

1 .0 

-15.2221 

1.0 

-12.9689 

S . 
J 

c . 

J 

formula 

H 

- 6.089 

H 

H2 

-17.164 

2H 

h2o 

-34.054 

2H  -t-  0 

N 

- 5-914 

N 

n2 

-24.721 

2N 

NH 

-14  . 9 36 

N + H 

NO 

-24.100 

N + 0 

0 

-10.708 

0 

°2 

-26.662 

20 

OH 

-22.179 

0 + H 

4.06  (io*c) 
1.48  ( 1.0" 1 ) 
7.63  ( io‘J ) 
l.4l  ( lO*3) 
4.83  ( io"J) 

-4 

6.93  ( '-0  ) 

2.72 

L.79  (iO'") 

3.73  (10'2' 
9.69  (10'^'j 


A. 4 SODA -POP  (Shapley  and  Cutler  [48]) 


Dimensions:  m = 9,  n = 17,  k = 2 

F*/RT  =■  -3128.86 


i 

B. 

1 

b. 

1 

* 
z . 
1 

1 

°2 

5.27583 

- 2.9421 

2 

co2 

6.07349 

-10.4964 

3 

N2 

82.58040 

-H.7II5 

4 

H+ 

52.81000 

-21.1660 

5 

oh" 

52.83950 

-18.2997 

6 

Cl" 

0.08005 

- 6.3746 

7 

+ 

Na 

0.08813 

- 6.2785 

8 

4 

K 

0.04829 

- 4.?834 

9 

GLUCOSE 

0.02000 

- 7.7625 

| 

i 

i 


A. 4 SCDA-POP  Cont. 


17 


GLUCOSE 


U.O 


formula 

t * 

Xi 

°2 

5-2798 

co2 

6.0427 

N2 

82 . 5800 

+ ~ 

II  + OH 

6.1075 

GAS  PHASE  x - 

100.006 

°2 

4.39  (jo"5) 

C°2 

1.50  (10‘5) 

N2 

3-85  (10‘4) 

H + OH 

40.  [0 

+ 

H 

3-02  do’8; 

oh" 

5.09  flO*4 

Cl’ 

8.00  ( 1.0" c i 

Na 

8.81  (1.0*3, 

4- 

K 

4.83  ( iC'c’j 

C02  + OH 

2.94  ( 10  J 

C02  +■  H+  + OH" 

1.83  ( 10”" , 

C0?  - H+  + OH’ 

4 . 86  ( I 0 ^ ) 

GLUCOSE 

2.00  ( iC'"’) 

LIQUID  PHASE  i 

46.97 

LIQUID  PHASE 


46.97 


A. 5 RESPIRATORY  SYSTEM  (Dantzig,  DeHaven  and  Sams  [17]) 


A. 5 RESPIRATORY  SYSTEM  Cont 


k i 

sj 

cj 

formula 

# 

XJ 

3 17 

°2 

0.0 

°2 

6.252  (10"5) 

' 18 

co2 

0.0 

C02 

6.888  (10"6) 

19 

N2 

0.0 

N2 

1.910  (10‘S 

20 

H+ 

0.0 

H+ 

1.276  (io"10) 

21 

oh“ 

0.0 

oh" 

3.929  (IO-5)  I 

22 

Cl" 

0.0 

Cl" 

8.140  (io"2) 

23 

K+ 

0.0 

K+ 

5.000  (10'2) 

+ 

24 

HgO 

-39.23 

H + OH 

23.29 

25 

HCO" 

3 

-21.20 

COg  + OH- 

1.858  (10"2) 

26 

¥*3 

0.0 

COg  + H+  + OH" 

6.273  (IO-25) 

27 

o°; 

6.25 

COg  - H+  + OH" 

4.095  (io-5)  ; 

28 

HBl" 

0.0 

HBl” 

2.951  (io“4) 

29 

HBlOg 

-16.23 

Og  + HBl" 

8.795  (IO-5) 

30 

HPr" 

0.0 

HPr" 

1.190  (io“2) 

RED  CELLS  = 

23.463 
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A. 6 PLASMA  MODEL  (jlasen  [14] ) 

Dimensions:  m = 1 6,  n = 21,  k = 1 

F*/RT  = -832.57 


Bi 

bi 

* 

z 

i 

R1  ; 

1.9073  (10-5) 

-13.9225 

R2 

1.1423  (10-2) 

-10.5336 

R3 

1.7458  /10"4) 

-11.7084 

R4 

21.108 

-21.9996 

R5 

21.119 

-18.5960 

r6 

4.2369  (10-2) 

- 6.2166 

R7 

5.8757  (io-2) 

- 5.8896 

R8 

1.8593  (io‘3) 

- 9.3^29 

R9  i 

7.3176  (10'4) 

-10.2754 

RIO 

3.9100  (io-4) 

-10.9021 

Rll 

1.1382  (io-4) 

-12.1362 

R12 

2.9362  (io-4) 

-II.6898 

R13 

7.1448  (lO~4) 

-10.2993 

Rl4 

1.1643  (10'5) 

- 9.8109 

R15 

1.0177  (IO-3) 

- 9.9455 

R16 

3.8600  (10"4) 

-10.9150 
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I 


A. 6 PLASMA  MODEL  Cont. 


3 

C3 

formula 

* 

X3 

1 

Cl 

0.0 

R1 

1.91  do"3) 

2 

C2 

0.0 

R2 

5.65  do'4) 

3 

C3 

0.0 

R3 

1.75  do"4) 

4 

C4 

0.0 

R4 

1.61  do"8) 

5 

C5 

0.0 

R5 

2.18  do"7) 

6 

c6 

0.0 

r6 

4.24  (10'2) 

7 

C7 

! o.o 

R7 

5.88  do"2) 

8 

C8 

0.0 

r8 

1.86  do"3) 

9 

C9 

0.0 

R9 

7.32  (10-4) 

10 

CIO 

0.0 

RIO 

3.91  (10'4) 

11 

Cll 

0.0 

Rll 

1.14  (lo'4) 

12 

C12 

0.0 

R12 

1.78  (10'4) 

13 

C13 

0.0 

R13 

7.14  (10-4) 

14 

Cl4 

0.0 

R14 

1.16  (10"3) 

15 

C15 

0.0 

R15 

1.02  (10"3) 

l6 

Cl6 

-21.35 

R2  + R5 

1.08  do"3) 

17 

C17 

-32.84 

R2  + R4  + R5 

8.04  do-7) 

18 

Cl8 

6.26 

R2  - R4  + R5 

1.46  <10"5) 

19 

C19 

-39-39 

r4  + R5 

21.1082 

20 

C20 

0.0 

Rl6 

3.86  do"4) 

21 

C21 

-20.57 

R4  + R12 

1.16  do"4) 

A. 7 FETUS  MODEL  (Clasen  [14] ) 


Dimensions:  m = 19,  n = 51>  k = 7 
f*/RT  = -1869.55 


i 

Bi 

bi 

1 

R1 

0.61550 

2 

R2 

0.28000 

3 

; R3 

37.370 

4 

R4 

46.115 

5 

R5 

0.00581 

6 

r6 

46.157 

7 

1 R7 

0.07925 

8 

R8 

0.07487 

9 

R9 

0.08292 

10 

RIO 

0.01022 

11 

Rll 

0.05720 

12 

R12 

0.00250 

13 

R13 

0.0 

14 

R12 

0.0 

; 15 

R15 

0.0 

16 

Rl6 

0.0 

17 

R17 

0.0 

18 

Rl8 

0.00295 

19 

R19 

0.01451 

- 2.092 
-10.610 
-11.792 
-22.005 

- 6.852 
-17.591 

- 4.219 
-59.995 

- 6.297 

- 6.242 

- 0.875 
-18.561 

2.171 

- 2.119 

- 5.846 

- 0.449 

- 0.200 
-11.827 
- 8.681 
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A. 7 FETUS  MODEL  Cont. 


formula 


Rl 

R2 

R3 

R4  + r6 


PHASEl 


R1 

R2 

R3 

R5  - R13 
R4  + R13 
r6  - RI3 
R7  - R13 
r8  + RI3 
R9  + R13 
r4  + r6 

BR  + R6  - R13 
R2  + R4  + r6 
R2  - R4  + r6  - 2R13 
RIO  - R13 
R13  + Rl8 
R13  + R1 


PHASE2  x2 


.6057 
.2646 
3.7367 
1 


4.9070 


1.43  fio" 


46.03 
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A. 7 FETUS  MODEL  Cont. 


formula 


21 

C21 

0.0 

R5 

22 

C22 

10.450 

R1 

23 

C23 

, 0.0 

R2 

24 

C24 

- 0.500 

R3 

25 

C25 

0.0 

R4 

26 

C26 

0.0 

r6 

27 

C27 

0.0 

R7 

28 

C28 

-39.389 

R8 

29 

C29 

0.0 

R9 

30 

C30 

3.373 

R4 

31 

C31 

-21.490 

R2 

32 

C32 

-32.840 

R2 

33 

C33 

6.120 

R2 

34 

C34 

0.0 

Rll 

35 

C35 

0.0 

R12 

36 

C36 

- 1.903 

R1  ■ 

37 

C37 

- 2.900 

2R1 

38 

C38 

- 3.362 

3R1 

39 

C39 

- 7.485 

4ri 

40 

C40 

2.606 

Rl8 

4l 

C4l 

0.0 

R19 

R2  + R4  + R6 
R2  - R4  + r6 


R1  + R12  - JRllf  - R15 
2R1  + R12  - 2R14  - 2R15 
3R1  + R12  - Rl4  - 3RI5 
4R1  + R12  - 4R15 


PHASE3 


1.48  do"4) 
4.90  (10-7) 
3.38  (10"6) 

1.71  (10“6) 

3.82  (io"11) 
3.84  <10-9) 
2.02  do"3) 

7.49  do"2) 
2.53  do"4) 

3.66  do"20) 
2.04  do"4) 

4.82  do"9) 
7.49  do"7) 

5.72  do"2) 

5.72  do’6) 

2.66  do"5) 

i 

5.01  do"5) 
5.51  do"5) 
2.36  do"3) 
7.40  (io-8) 
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formula 


R4  + Rl4 
Rl4 

-R4  + Rl4  - 4ri6 
PHASE4  5 

R4  + R15 
R15 

■R4  + R15  - 4R17 
PHASE5  xc  = 

Rl6 

Rg  - r4  + Rl6 
PHASE6  : 


A. 8 SEA  POWER  (Duffin,  Peterson  and  Zener  [25,  p.  127]) 
Dimensions:  m = 7>  n = 10,  k = 4 


k 

J 

Exponents 

* 

uj 

0 

1 

1.0 

A 

126.7476 

OBJECTIVE  = 

126.7476 

1 

2 

40.000 

^v^SHHSSBSI 

1.0000 

*i- 

0.8265 

2 

3 

0.1800 

Rmn 

1.000 

"2  - 

0.4205 

3 

4 

44.5000 

Q"1  * Q_1 

0.7753 

5 

6.00  (lo‘8) 

5 -1  -1 

A*tr*Q  *a 

0.1915 

6 

2.15  do"8) 

u2  * cf1  * r"1 

0.0332 

.....  _V 

1.2898  ; 

4 

7 

1.0 

a 

0.5000 

8 

1.0 

P 

0.3204 

9 

1.0 

P’ 

0.1630 

10 

1.0 

r 

0.0166 

2.5795 
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A. 8 SEA  PCWER  Cont. 


Symbol 


Low  bd. 

Upper  bd. 

1.000 

500.0 

126.7476 

1.000 

300.0 

114.7878 

10.000 

1000.0 

113.0615 

0.001 

1.0 

0.5000 

0.001 

1.0 

0.3204 

0.001 

i.o  ; 

0.1630 

0.001 

1.0 

0.0166 

A. 10  CONDENSER  DESIGN  (Avriel  and  Wilde  [2]) 
Dimensions:  m-4,  n=8,  k=2 


Exponents 


•GO 

DO*1  * L* ^ * N*®/7 

98.6717 

'90 

mV5.  L-l  » „-l/5 

171.0548 

1.570 

DO  * L * N 

410.6575 

0.0382 

-kT"^»8  u * jt  «t*  1*8 

DI  * L * N 

54.7279 

i80 

Dl" 1 * L*1  * N*1 

162.6811 

OBJECTIVE  = 

897.8829 

(10-3) 

DO*1 

0.0980 

-1 

DI  * DO 

0.9020 

>- 

i-j 

u 

0.3563 

DO 

1.0000 

II 

0.0089 

Symbol 

Low  bd. 

Upper  bd. 

* 

DI 

0.02 

0.080 

0.0752 

DO 

0.03 

0.084 

O.O833 

L 

1.00 

200.0 

28.0715 

N 

1.00 

1000.0 

Hl.8138 

NOTE:  Dl  and  DO  are  each  considered  a single  symbol. 
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